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CHAPTER  I 
INTRODUCTION 


This  study  describes  the  performance  of  a  Kalman  filter  imple¬ 
mentation  of  a  sample  data  delay  lock  loop.  The  (non-sampled)  delay 
lock  looo  was  originally  studied  by  Spi lker[l 7,18]  and  Gill[16]  for 
use  in  the  tracking  of  pseudo-noise  (PN)  coded  signals.  A  sampled 
data  version  of  this  loop  was  developed  by  Huff  et  al . [19,20,21]  for 
receiver  timing  in  a  time  division  multiple  access  (TDMA)  satellite 
communications  system. 

To  facilitate  receiver  timing  in  this  system,  a  PN  coded  clock 
pulse  is  periodically  transmitted  from  the  satellite.  In  the  receiver 
the  sample  data  delay  lock  loop  "tracks"  the  clock  pulse  by  control¬ 
ling  the  frequency  and  phase  of  a  local  square  wave  oscillator  (local 
clock).  In  the  interval  between  clock  pulses,  the  transitions  of 
the  local  clock  are  used  to  time  the  correlators  in  the  receiver  de¬ 
tection  circuits.  Errors  in  clock  tracking  (biases  and/or  jitter) 
result  in  imperfect  timing  o 
crease  in  bit  error  probabil 
mize  the  bias  and  error  vari 
desire  to  increase  PN  code  ratios  to  40  MHz  makes  the  timing  require¬ 
ment  more  critical.  Furthermore,  the  incorporation  of  highly  maneu¬ 
verable  terminals  ( i . e . ,  aircraft)  in  the  system  compounds  the  track¬ 
ing  problem  by  adding  time  variable  doppler  effects. 

In  order  to  achieve  the  desired  tracking  accuracy  the  Kalman 
filter  algorithm  was  chosen[6].  This  algorithm  is  known  to  yield 
both  oDtimum  estimates  (unbiased,  minimum  variance  estimates)  and 
a  measure  of  its  own  performance,  an  error  covariance  matrixll,?,3,4] 
In  order  to  be  implemented,  however,  the  Kalman  filter  does  reauire 
a  mathematical  model.  This  is  a  set  of  equations  which  describes 
the  process  to  be  filtered.  The  error  variances  contained  in  the 
error  covariance  matrix  are  the  error  variances  for  the  modeled  (as 
opnosed  to  the  actual)  situation.  One  should  have  some  knowledge 
concerning  model  accuracy  before  basing  a  design  upon  the  error  vari¬ 
ances  of  the  model.  The  model  chosen  for  the  tracking  filter  appli¬ 
cation  is  a  statistical  representation  of  the  time  variable  propagati 
delays  generated  by  terminal  motion,  specifically  aircraft  flight. 
Checks  on  the  accuracy  of  the  model  would  involve  comparisons  of  the 
error  variances  yielded  by  the  algorithm  with  the  error  variances 
actually  obtained  by  implementing  the  filter  in  a  moving  terminal. 
This  is  beyond  the  scope  of  this  study. 


f  the  detectors  with  a  corresponding  in- 
ityl.7,22].  Thus,  it  is  desirable  to  mini 
ance  (jitter)  of  the  local  clock.  The 
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The  purpose  of  this  study  is  to  evaluate  the  performance  of  a 
Kalman  tracking  filter  independently  of  the  question  of  model  accuracy. 
Given  the  mathematical  model,  the  structure  and  gains  of  the  filter 
are  determined  by  the  algorithm.  The  resulting  filter  may  then  be 
analyzed  for  its  response  to  deterministic  inputs  and  noise  using 
conventional  techniques.  Thus,  the  model  is  used  only  to  determine 
the  filter  structure  and  not  to  determine  performance.  Because  the 
Kalman  filter  is  in  state  variable  form,  the  most  convenient  techni¬ 
que  for  analysis  is  computer  simulation.  Other  techniques  include 
state  variable  techniques  and  z-transform  methods.  Since  the  Kalman 
filter  is  linear,  suoerposi tion  ma”  be  used  to  determine  the  response 
to  linear  combinations  of  inputs. 

ChaDter  II  of  this  study  presents  an  overview  of  the  SDDLL,  Kalman 
theory,  the  model  chosen  for  this  application,  and  some  discussion 
of  computer  simulation.  In  Chapter  III  the  steady-state  response 
of  the  Kalman  filter  to  input  noise  is  analyzed  using  the  state  vari¬ 
able  techniques.  Monte  Cairo  simulation  is  used  to  verify  the  results. 
In  addition,  an  analysis  using  z-transform  techniques  is  presented. 

The  steady-state  response  of  the  Kalman  tracking  filter  to  determin¬ 
istic  delay  changes  caused  by  a  terminal  maneuvering  through  180° 
turns  is  presented  in  Chapter  IV.  In  Chapter  V  the  combined  response 
to  noise  and  maneuvers  is  considered.  Three  examples  are  presented 
which  illustrate  the  usefulness  of  the  preceeding  data.  Chapter  VI 
documents  an  analysis  of  the  lockup  characteristics  of  the  tracking 
loop.  Here  the  response  to  step,  ramp,  and  parabolic  delay  changes 
is  determined  for  Kalman  tracking  filters  in  both  the  transient  and 
the  steady  state.  The  number  of  iterations  required  for  lock-up  (con¬ 
vergence)  for  special  cases  of  interest  is  also  determined.  Chapter 
VII  presents  the  conclusions  of  the  study. 
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CHAPTER  II 

KALMAN  FILTERING  AND  THE  TOMA  CLOCK  LOOP 


A.  Introduction 


In  this  chapter  a  short  review  of  the  TDMA  clock  loop  is  presented. 
This  is  followed  by  an  overview  of  Kalman  filter  theory  for  the  pur¬ 
pose  of  introducing  the  notation,  equations,  and  concepts  to  be  used 
in  subsequent  analyses. 

B.  The  TDMA  Clock  Loop 

The  OSU/RADC  TDMA  system  uses  pseudo-noise  (PN)  coding  for  im¬ 
proved  multipath  and  anti-jam  performance!.^] .  A  pulsed  envelope 
PN  coded  "clock"  pulse  is  periodically  transmitted  on  the  down-link 
to  facilitate  PN  code  tracking  at  the  receiver.  Figure  1  shows  the 
envelope  of  the  clock  pulse.  Here,  A  is  the  chip1  duration,  M  is 
the  number  of  chips  per  clock  pulse,  T  is  the  duration  of  the  clock 
pulse,  and  T^  is  interval  between  leading  edges  of  successive  clock 
pulses. 


Figure  1.  The  envelope  of  the  ospudo-randcm  network  clock  signal. 


^A  chip  is  the  elementary  unit  of  a  PN  code  and  is  analogous  to  a 
hit  in  a  binary  data  stream. 
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Tracking  of  the  PN  code  is  accomplished  at  the  receiver  by  a 
sampled-data  delay-lock  loop  fSDDLL)[16, 17,18,19,20,21] ,  a  block  dia¬ 
gram  of  which  is  shown  in  Figure  2.  The  function  of  this  tracking 
loop  is  to  phase  lock  a  locally  generated  replica  H.e.,  the  "local 
clock")  of  the  clock  pulse  with  the  down  link  signal.  The  local  clock 
siqnal  is  generated  by  a  pulse  generator  and  feedback  shift  register. 
Delayed  and  advanced  versions  of  the  local  clock  are  correlated  with 
the  down-link  clock  siqnal  in  the  upper  and  lower  mixer- bandpass  fi lter- 
detector  circuits,  respectively.  The  difference  of  the  correlator 
outputs  is  then  sampled,  yielding  the  sampled  error  voltage  E  .  This 
is  shown  in  Figure  3  as  a  function  of  the  timing  error  (e)  between 
the  received  and  local  clock  pulses.  The  error  voltage  is  then  fil¬ 
tered  by  the  discrete  loop  filter.  This  filter  can  be  implemented 
with  a  variety  of  algorithms.  The  present  system  uses  a  simple  av¬ 
eraging  filter.  This  study  considers  implementation  with  a  Kalman 
algorithm.  The  output  of  the  discrete  filter  is  then  fed  to  the  clock 
correction  circuitry  to  adjust  the  time  base  of  the  local  clock. 


Kalman  Filter  Theory 


The  Kalman  estimator  algorithms  are  optimum  algorithms  in  that 
they  yield  minimum  variance,  unbiased  estimates  for  gaussian  processes. 
Unlike  the  Wiener  filter,  which  produces  optimum  estimates  in  the 
steady-state,  the  Kalman  algorithms  produce  optimum  transient  as  well 
as  optimum  steady-state  response.  They  also  yield  optimum  estimates 
for  non-stationary  or  periodically  stationary  processes  for  which 
a  steady-state  does  not  exist.  Because  of  the  feature  of  optimality 
the  Kalman  algorithms  were  chosen  as  a  method  for  potentially  upgrading 
the  performance  of  the  SDDLL. 


These  algorithms  are  cast  in  state  variable  form  and  require 
mathematical  models  of  both  the  random  process  to  be  filtered  and  of 
the  observations  of  the  process.  The  random  process  should  be  modeled 
in  the  form 


z(k+i )  =  $(k)z(k)  +F(k)u(k)  (D 

where  k  is  the  iteration  number  and  where 
z(k)  =  n  dimensional  state  vector, 

u^k)  =  r  dimensional  zero  mean  white  process  noise  vector, 

T(k)  =  nxr  dimensional  forcing  matrix, 
and  ^ 

4>(k)  =  nxn  dimensional  state  transition  matrix  operating  on  the 
kth  state 
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FEEDBACK  LOGIC  NOT  SHOWN 


Figure  2,  Block  diagram  of  the  sampled  data  delay  lock  loop  (SDDLL) 


Figure  3.  Sampled  error  voltage  (Es)  as  a  function  of  timing  error(c) 
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The  observation  model  should  be  in  the  form 


x ( k )  =  H(k)  z(k)  +  v(k) 


(2) 


where 

and 


x(k)  j  s  dimensional  observation  vector, 
vfk)  =  s  dimensional  zero  mean  White  measurement,  noise 
vector, 

H f k )  -  sxn  dimensional  observation  matrix. 


Together  these  two  equations  model  the  input  to  the  algorithms.  The 
first  equation  qenerates  the  state  vector  z.  From  Equation  (21  the 
observation  of  the  process  is  composed  of  a  linear  combination  of  the 
states  (H(k)z(kU  plus  additive  white  measurement  noise  (vfk)).  The 
covariance  matrices  of  the  two  noise  vectors  and  v(k)  are  also  re¬ 
quired  for  implementation  of  the  algorithm.  These  are  defined  by 

E{u(j)uT(k)}  =  V  (k),  the  rxr  dimensional  process  noise 

covariance  matrix, 

and  T  A 

E { v ( j ) u  ( k ) }  =  V  (k),  the  sxs  dimensional  observation 

noise  covariance  matrix. 


Here  E  is  the  expectation  operation  and  T  represents  transposition. 
Thus,  the  models  required  for  implementation  of  the  Kalman  algorithms 
are  completely  determined  by  specifying  the  matrices  $(k),  r(k),  H(k), 
Vu(k),  and  Vy(k ) . 

The  two  Kalman  algorithms  investigated  for  use  in  the  SDDLL  are 
the  Kalman  one  step  predictor  algorithm  and  the  Kalman  filter  algo¬ 
rithm.  The  one  step  predictor  equation  is  given  byU.,2,3,4] 

z(k+l/k)  =  <f(k)2(k/k-l)+Kp(kl[x(k)-H(klz(k/k-l)]  ,  (3) 


where 


z(k+1/k)  -  the  n  dimensional  prediction  of  the  state  vector  at 
the  (k+l)st  iteration  given  k  data  samples 


Kp(k)  =  the  nxs  dimensional  predictor  gain  matrix. 


This  algorithm  is  iterative;  in  order  to  generate  the  prediction  of 
the  (k+l)st  state,  the  prediction  of  the  kth  state  is  required.  Thus, 
to  start  the  algorithm,  an  initial  estimate  must  be  obtained  by  means 
other  than^Equation  (3).  Also,  notice  that  the  term  in  brackets 
(x(k)-H(k)z(k/k-l))  is  the  observation  minus  the  prediction  obser¬ 
vation.  This  term  is  sometimes  called  the  innovation  and  represents 
the  new  information  contained  in  the  measurement.  In  the  present 
application  the  observation  will  be  modeled  as  noise  corrupted  propa¬ 
gation  delay  from  a  satellite.  A  measurement  of  the  propagation  delay 


is  impossible  in  the  SDDLL,  which  can  only  measure  the  difference 
in  the  time  bases  of  the  downlink  and  local  clock  signals.  However, 
bv  forcing  the  local  clock  to  follow  the  predicted  delay,  a  measure¬ 
ment  of  the  term  in  brackets  (i.e.,  delay  minus  predicted  delay,  the 
innovation!  is  obtained1.  Equation  (3)  states  that  the  new  pre¬ 
diction  (z(k+1/k))  is  obtained  by  advancinq  the  old  prediction 
fz(k/kl) )  one  step  ahead  with  the  state  transition  matrix  ($(k))  and 
addinq  weighted  (by  the  gain  matrix,  K  (k)  amounts  of  new  information 
(innovation).  p 

The  gain  matrix  needed  for  use  in  Equation  (3)  is  given  by 

K  (k)  =  3>(k)P(k)HT(k)[H(k)P(k)HT(k)+Vv(k)]_1  ,  (4) 

Here,  <5(k),  H(k)_,  and  Vy(k)  were  specified  for  the  process  and  obser¬ 
vation  models.  P(k)  isvthe  prediction  error  covariance  matrix,  de¬ 
fined  by 

P(k)  -  E{z  (k)zT(k)},  the  nxn  dimensional  prediction  error 
e  covariance  matrix, 

where 

z  (k)  =  z(k)-z(k/k-l),  the  one  step  prediction  error  at  the  kth 
e  point  given  k-1  data  samples. 

The  error  matrix  may  be  calculated  iteratively  via  the  equation 

P(k+1 )  =  4>(k)P(k)4>T(k)+r(k)Vu(k)rT(k)  (5) 

-  4>(k)P(k)HT(k)[H(k)P(k)HT(k)+Vv(k)]'1  H  ( k )  P  ( k )  ^T(  k ) . 

This  equation  is  a  function  of  $(k),T(k),  H(k),  V  (k )  and  V  (k),  all 
of  which  are  specified  for  the  process  and  observation  models.  Since 
Equation  (5)  is  iterative,  an  initial  error  covariance  matrix  is  needed 
for  the  first  iteratiofi.  This  matrix  may  be  obtained  from  an  error 
analysis  of  the  initial  estimate  used  to  begin  the  iterations  of  Eq¬ 
uation  (3).  If  the  process  and  observation  models  are  stationary 
(non  time  varying)  the  error  covariance  matrix  (and  hence  also  the 
predictor  gain  matrix)  will  approach  asymptotic  or  steady-state  values 
as  the  iterations  proceed.  In  addition  to  its  use  for  computing  the 
gain  matrix  (Equation  (4)),  the  error  covariance  matrix  also  provides 
a  measure  of  the  predictor  performance  in  response  to  the  modeled 
process.  The  diagonal  elements  of  this  matrix  are  the  variances  of 

^To  avoid  confusion,  the  term  observation  will  be  used  exclusively 
to  represent  x(k)  and  the  term  measurement  will  be  used  exclusively 
to  represent  the  error  signal  obtained  from  the  SDDLL  (i.e.,  the  in¬ 
novation). 


the  prediction  error  vector  z£.  If  the  process  model  accurately  re¬ 
flects  the  physical  situation,  then  the  error  convariances  obtained 
in  an  actual  implementation  will  approach  those  given  theoretically 
bv  Equation  (5). 

Equations  (3),  14),  and  15)  are  collectively  known  as  the  Kalman 
one  step  predictor  algorithm.  Equation  (3)  forms  the  actual  prediction. 
A  block  diagram  of  the  equation  is  shown  in  Figure  4.  Equations  (4) 
and  (5)  are  used  for  the  calculation  of  the  gain  matrix  used  in  Eq¬ 
uation  (3).  These  equations  can  be  used  in  either  of  two  methods. 

Figure  5  shows  a  block  diagram  for  implementing  the  entire  prediction 
algorithm  in  real  time.  After  updating  the  prediction,  the  error 
covariance  and  gain  matrices  are  calculated  prior  to  taking  the  next 
measurement.  The  measurement  is  then  obtained,  the  prediction  up¬ 
dated,  and  the  process  repeated.  This  method  is  possible  as  long 
as  the  time  interval  between  measurements  is  long  enough  to  allow 
the  necessary  computations.  Alternately,  since  the  error  covariance 
and  gain  matrices  are  functions  only  of  the  process  and  measurement 
models  and  not  of  the  actual  measurements,  they  may  both  be  computed 
and  stored  prior  to  enabling  the  predictor  equation.  With  the  gain 
matrix  for  each  iteration  available  in  memory,  a  considerable  savings 
in  computation  time  can  be  achieved  for  the  operation  of  the  predictor 
equation. 

The  optimal  filtered  estimate  is  related  to  the  optimal  one  step 
prediction  bv[l ,?] 

zfk-1)  =  $-1 (k)z(k/k-l )  16) 

or 

z(k/k-l)  =  <Kk-l)z(k-l)  (7) 

where  $_1(k)  may  be  thought  of  as  the  reverse  state  transition  matrix 
(i.e.,  the  transition  matrix  going  from  the  kth  state  to  the  (k-l)st 
state)  and  where 

A  A 

z(k-l)  -  the  (n  dimensional  vector)  filtered  estimate  of 
the  state  z(k-l)  after  processing  the  (k-1)  data 
sample. 

Equations  (3),  (6),  and  (7)  may  be  used  to  derive  Equation  (8),  the 
Kalman  filter  equation: 

z(k)  =  <J>(k-l)z(k-l)+Kf(k)Lx(k)-H(k)$(k-l)z(k-l)].  (8) 

Here  Kf(k)  is  the  filter  gain  matrix  for  the  kth  data  point.  Note 
that  tne  term  $(k-l)z(k-l)  is  the  one  step  prediction  z (k/k-1 ) .  Fur¬ 
thermore,  the  term  inside  the  brackets  is  again  the  innovation.  Eq¬ 
uation  (8)  states  that  the  filtered  estimate  is  obtained  from  a  sum 
of  two  components.  The  first  component  is  the  previous  filtered  es¬ 
timate  advanced  one  step  forward  by  the  state  transition  matrix.  The 


z(k+l/k  )  =  <1>  z  (  k/k-l )  +  Kp(k)  [x(k)~  Hz*(  k/k-l)] 

Figure  4.  Block  diagram  of  the  Kalman  Predictor  Equation 
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second  term  is  the  innovation  (new  information)  weighted  by  the  filter 
adin  matrix  K^(kl . 

The  filter  gain  matrix  is  qiven  by 

Kf(k)  =  <t>_1(k)Kp(k)  =  P(k)HT(k)LH(k)P(k)HT(k)+Vv(k)]'1'.  (9) 

As  in  the  case  of  the  predictor  gain  matrix,  the  filter  gain  matrix 
is  also  a  function  of  the  prediction  error  covariance  matrix.  The 
prediction  error  covariance  matrix  is  calculated  iteratively  in  Eq¬ 
uation  (5),  Equations  (8),  (9),  and  (5)  collectively  from  the  Kalman 
filter  algorithm[l,2] .  These  equations  are  used  in  a  manner  analo¬ 
gous  to  the  predictor  algorithm  with  Equations  (5)  and  (9)  being  used 
to  compute  the  gain  matrix  necessary  for  the  computation  of  the  fil¬ 
tered  estimate  (Equation  (8)). 

The  error  covariance  matrix  for  the  filtered  (as  opposed  to  the 
predicted)  estimate  is  defined  by 

P(k)  =  E[z  (k)S(k)] 

where 

z£(k)  =  z(k)-z(k) 

is  the  n  dimensional  filter  error  vector  at  the  kth  iteration.  The 
filter  error  covariance  matrix  may  be  obtained  from  the  predictor 
error  covariance  matrix  via  the  equation 

P(k)  «  [l-Kf(k)H(k)]p(k)  .  (10) 

The  matrix  P(k)  is  a  measure  of  filter  performance  just  as  the  matrix 
P(k)  is  a  measure  of  predictor  performance.  Both  P(k)  and  P(k)  are 
error  covariance  matrices  for  the  modeled  situation.  The  values  given 
by  these  matrices  may  be  considered  accurate  only  if  the  model  is 
accurate  in  its  representation  of  the  actual  process.  Thus,  the  ac¬ 
curacy  of  a  given  process  model  must  be  established  before  the  error 
covariance  matrices  can  be  considered  reliable. 

Fortunately,  it  is  possible  to  evaluate  the  performance  of  the 
Kalman  algorithms  without  establishing  the  accuracy  of  the  process 
model.  Suppose  the  process  and  observation  models  are  specified. 

This  determines  the  structure  and  gain  of  the  resulting  predictor 
and  filter  equations  (Equations  (3)  and  (8),  respectively).  The  re¬ 
sponses  of  these  equations  to  various  inputs  (such  as  steps,  ramps, 
white  noise,  etc.)  can  then  be  determined  using  conventional  tech¬ 
niques.  Note  that  the  responses  obtained  can  be  considered  true  re¬ 
sponses  even  if  the  process  model  does  not  accurately  represent  the 
random  process  (i.e.,  even  if  there  are  modeling  errors).  This  is 
so  because  the  structure  and  gain  of  the  Kalman  filter  and/or  pre¬ 
dictor  are  functions  only  of  the  model,  not  of  the  model  accuracy. 
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for  example,  suppose  a  process  model  is  chosen  which  accurately 
models  the  propagation  delays  generated  by  a  slowly  moving  mobile 
terminal  I  such  as  an  automobile).  Using  this  process  model,  the  step 
response  of  Equation  (3)  I  the  predictor  equation)  is  determined. 

Now,  suppose  the  same  process  model  is  used  to  represent  the  propa¬ 
gation  delays  generated  by  a  rapidly  moving  mobile  terminal  (such 
as  an  aircraft).  Obviously,  some  modeling  errors  will  probably  be 
present.  However,  the  step  response  will  be  the  same  simply  because 
the  same  process  model  is  used. 

This  technique  is  used  in  the  remainder  of  this  study  to  obtain 
performance  data  for  a  wide  variety  of  inputs.  All  of  the  data  pre¬ 
sented  in  the  sequel  can  be  considered  valid  regardless  of  the  process 
and  observation  model  accuracies.  This  is  generally  not  true  of  the 
performance  data  contained  in  the  previous  study fc  J,  most  of  which 
is  valid  only  if  the  process  and  observation  models  are  accurate. 

0.  The  Process  and  Observation  Models 
for  the  Clock  Loop 

The  Drocess  model  for  the  Kalman  clock  loop  filter  is  chosen 
from  Sinqer  and  Behnke[5]  and  models  delay  and  the  first  two  deriv¬ 
atives  of  delay  with  respect  to  time  as  a  random  process.  Specif¬ 
ically, 


Vk+1)  =  Vk)+f0(k)VVk)Tf/2’ 


(Ha) 


f0(k+l)  -t0(k)+ToTf, 


(lib) 


T0(k+1)  =  PTo(k)  +  u(k) 


(11c) 


where 


T  (k)  =  =  downlink  propagation  delay  at  the  kth  sample 

'  c  instant, 

A  dVk> 

f(k)  -  — jt -  =  downlink  propagation  delay  rate  at  the  kth 

1  '  sample  sample  instant  (delay  velocity), 

A  d?Vk) 

TOO  -  - p —  =  downlink  propagation  delay  acceleration  at 

J  dt^  the  kth  sample  instant. 


c  -  speed  of  light, 

r(k)  -  downlink  range  at  kth  sample  instant, 


and 


P  =  E {T ( t+r )  •  T( t ) } j  _T  =  delay  acceleration  correlation 

T"  f  coefficient. 

Equations  (11a)  and  (lib)  are  Taylor's  series  representations  for 
delay  (three  term  Taylor  series)  and  delay  rate  (two  term  Taylor 
series).  Equation  (11c)  models  acceleration  as  a  correlated  random 
process.  The  coefficient  is  chosen  from  an  exponential  correlation 
function  so  that 


p(  t) 


_  e~At  = 


=  1-Ax  + 


<AiY 


(12) 


where  A  has  units  of  sec"^  and  is  considered  to  be  the  inverse  of 
the  averaqe  maneuver  duration[5,f] .  For  t=T^  and  AT^  <<  1,  p  is  given 
approximately  by 

p=  1  -  ATf  .  (13) 

Equations  (11a),  (lib),  and  flic)  are  in  the  form  of  Equation  (1) 
where 


z(k)  = 


( \w\ 

Vk> 

y°(kv 

( i 


4>fk)  =  4>  = 


V 


=  state  vector. 


(14a) 


0  1 

n  0 


T^/2 


\ 


J 


=  constant  state  transition 
matrix,  fT4b) 


/ 


rfk)  =  r  = 


o 


0  J  =  constant  forcing  matrix, 

\  1 


M4c) 


12 


and 


u(k)  =  scalar  white  process  noise  with  constant  variance  a2 

(14d) 

Note  that  the  vector  u(k)  in  Equation  (1)  is  of  dimension  one  here 
(i.e.,  a  scalar).  Hence,  the  covariance  matrix  V  (k)  in  the  Kalman 
algorithms  (Equations  (5)  and  9))  reduces  to  the  Scalar  cr.  Squaring 
both  sides  of  Equation  (11c)  and  taking  expected  values  yields 

o|(  1-P?)  =  o2  .  (15) 

Thus,  au  is  easily  found  given  p  and  c^.  To  find  at  the  probability 

density  function  shown  in  Figure  6  is  assumed[5,6] .  The  quantity 
At  in  this  figure  represents  the  maximum  delay  acceleration  given 
by 

At  =  A/c  (16) 

where 

o 

A  =  maximum  terminal  acceleration  (m/sec  ), 
and 

c  =  speed  of  light  (m/sec). 

cr2  =  A*  (  I +  4P  -  P  )  /  3 
j  T  mo 


a 

P°  1/ 2  A  T 

Figure  6.  The  Probability  Density  Function  Function  for 
Dp  lay  Acceleration 
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From  this  density  function,  o^.-  is  easily  found  to  be 


?  ^ T 

cl  =  imp -p  i 
id  mo 


H7) 


where 


Pm  =  probability  of  maximum  acceleration 
and 

PQ  =  probability  of  zero  acceleration. 

The  observation  is  modeled  as 

x(k)  =  T0(k)  +  v(k).  (18) 

This  is  in  the  form  of  Equation  (2)  where  x(k)  is  now  a  scalar  and 
where 


v(k)  =  scalar  white  measurement  noise  with  stationary  variance 


(19a) 


and 


Hlk)  =  H  =  (1  0  01  =  constant  observation  matrix. 


(19b) 


is  a  scalar,  the  covariance  matrix  V  (k)  in  Equations  (4), 
i  18)  becomes  cr.  In  Reference  7,  Equation  1391  it  is  showi 


Since  v 

IS),  and  18)  becomes  In  Reference  7,  Equation  1391  it  is  shown 
that  the  variance  of  asingle  unprocessed  measurement  (samplel  in 
the  SODLL  is 


a?  .  1 


where 


1 .21 6  .  „  1.476 


(20) 


E  -  received  energy  in  the  clock  pulse  1 joules), 


and 


Nq  -  one  sided  noise  power  spectral  density  (watts/Hzl, 


C  =  < 


1  for  linear  detection 
[2  for  square  law  detection. 


Thus,  can  be  obtained  given  a  knowledge  of  E/NQ. 
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Equations  Cl),  (21,  fUa.b.c.d),  and  (19a, b)  together  with  vari¬ 
ances  (?  and  <?  completely  define  the  model  for  the  Kalman  filter 
considered  subsequently.  This  model  is  a  function  of  four  parameters: 

Tf,  <?  <?  and  X  for  p) ,  Neither  the  equations  nor  the  parameters 

vary  withtime;  hence  the  model  is  stationary.  A  typical  aircraft 
flight  might  consist  predominantly  of  mild  maneuvers  with  perhaps 
more  severe  manuevers  occasionally.  In  other  words,  typical  aircraft 
flight  is  non-stationary.  This  presents  problems  to  the  designer. 

If  the  selection  of  model  parameters  is  based  upon  average  maneuvers, 
the  tracking  error  may  become  excessively  large  during  severe  maneuvers. 

Alternately,  X  and  o^-  (and  hence  c?)  may  be  selected  to  be  representative 
of  more  severe  maneuvers,, perhaps  with  some  performance  degradation 
for  less  severe  maneuvers1.  Furthermore,  even  if  a  maneuver  is  speci- 
cified,  an  appropriate  selection  of  the  statistical  parameters  X  and 
2 

cTj  to  optimally  represent  the  maneuver  is  not  obvious.  To  confidently 

select  usable  model  parameters,  the  designer  needs  to  know  either 

(1)  the  relationships  of  the  model  parameters  to  the  performance  of 
the  aircraft  or  (2)  the  relationships  of  the  model  parameters  to  the 
filter  behavior  when  subjected  to  a  wide  variety  of  inputs. 

E.  Computer  Simulation 

The  problem  of  simulating  a  filter's  response  to  an  input  signal 
may  be  broken  into  two  parts,  fl)  simulating  the  input  signal  and 

(2)  simulating  the  filter.  In  the  Kalman  filter  SDDLL  described  in 
Reference  6  the  time  base  of  the  local  clock  serves  the  function 

of  both  the  optimal  estimate  and  the  optimal  prediction  of  the  delay 
function  T  (t).  Immediately  after  the  processing  of  a  new  sample, 
the  time  base  of  local  clock  serves  as  the  filtered  estimate.  Be¬ 
tween  data  samples,  the  time  base  is  the  optimum  prediction.  Immed¬ 
iately  before  the  next  data  sample  the  local  time  base  serves  as  the 
optimal  one-step  prediction.  This  is  the  "worst  case"  point  for  the 
filter,  the  time  at  which  the  largest  tracking  errors  generally  occur. 

To  investigate  the  worst  case  behavior  of  the  filter,  the  computer 
programs  used  in  the  remainder  of  this  report  simulate  a  Kalman  one- 
step  predictor.  This  is  easily  accomplished  by  direct  simulation 
of  Equations  (3),  (4),  and  (5). 

Simulation  of  the  predictor  during  transient  studies  requires 
an  initial  prediction  vector  z(0/-l)  and  an  initial  covariance  matrix 
P(o). 


It  should  be  noted  that  an  adaptive  scheme,  perhaps  with  accelerom¬ 
eter  inputs,  could  be  used  to  adjust  the  model  parameters  to  maintain 
an  optimal  filter  for  a  wide  variety  of  maneuvers.  This,  however, 
is  beyond  the  scope  of  this  report. 
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For  steady-state  studies,  the  covariance  matrix  may  be  iterated 
using  Equation  (5)  until  steady-state  is  achieved.  The  predictor 
qain  vector  may  then  be  computed  using  Equation  (4).  Alternately, 
one  may  compute  the  steady  state  gain  vector  and  covariance  matrix 
using  the  equations  given  in  Reference  6.  The  later  was  not  done  be¬ 
cause  (1)  the  computer  programs  described  subsequently  were  intended 
to  be  general  ( i . e. ,  used  both  for  transient  and  steady-state  simu¬ 
lations)  and  (2)  solution  by  iteration  serves  as  a  useful  check  on 
the  analytic  solutions  given  in  Reference  6. 


CHAPTER  III 

STEADY-STATE  RESPONSE  TO  NOISE 


Since  the  Kalman  predictor  is  linear,  superposition  may  be  used 
to  find  the  response  to  a  sum  of  a  deterministic  input  plus  noise. 
That  is,  one  may  find  the  response  due  to  noise  alone  and  then  sep¬ 
arately  find  the  response  due  to  the  deterministic  input  (without 
noise).  The  total  response  is  the  sum  of  the  individual  response. 
This  technique  was  used  in  this  study,  and  this  chapter  documents 
the  analysis  of  the  response  to  noise. 

A.  State-Variable  Approach 

In  state-variable  terminology,  the  solution  of  the  predictor 
response  to  measurement  noise  alone  is  accomplished  by  constraining 
the  state  vector  to  be  zero  for  all  iterations.  The  prediction  error 
is  then  given  by 

7  fk)  =  -  z  (k/k-1)  (21) 

n  n 

where  the  subscript  n  indicates  the  contribution  due  to  noise  only. 
Thus,  under  the  preceding  constraint,  the  error  covariance  is  the 
prediction  covariance.  That  is 

E \i  (k)  •  f  (k)}  =  E{?n(k/k-1)zjfk/k-l)}  (22) 

n  n 

The  prediction  is  given  by  the  recursion  relation  (from  Equation  (3)) 
z(k+l/k)  =  ($-k  (k)H)  z(k/k+l)+Kp(k)x(k)  (23) 

where  x(k)  is  the  observation. 

For  the  purpose  of  computing  the  gains,  the  problem  is  modeled 
normally  by?specifing  the  matrices  given  in  Equation  (14)  and  the 
variances  cr  and  cr.  When  computing  the  prediction,  the?observation 
x(k)  is  composed  only  of  white  noise  (E)  with  variance  of.  Here  we 
are  allowing  for„the  possibility  that  the  noise  may  not  have  the  mod¬ 
eled  variance  (of,)  but  instead  may  have  a  different  variance  (of). 

€,  and  v  both  represent  measurement  noise,  but  £  is  necessary  as^a 
dummy  variable  to  allow  for  noise  variance  other  than  those  modeled. 
Equation  (231  then  becomes 


Zn(k+l/k)  =  (*-k  (k>H)zn(k/k+l)  +  K  (k)t;(k)  (24) 
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Under  the  preceding  assumptions,  the  error  covariance  matrix  fi.e., 
F  ff  fk+1)>zr  (k+lil)  is  qiven  bv 

Pnrk+l)=(<I>-Kp(k)HiPn(k)($-Kp^iN)T  +  Kpo^Kp(k)  . 

When  E  is  stationary,  P  (k+1)  equals  P  (k)  for  sufficiently  large 
k.  Thus,  in  the  steadynstate 

P„  •  (j-KpH)Pn(^-KpH)T  *  0eVJ  (25) 

where  P  =  lim  P  fk)  and  K  =  lim  K(k). 

n  k-voo  n  p  k+*>  p 

Let  the  elements  of  P  be  qiven  by 

n  n  j 


">tn 

P?n  P3n 

V 

P2n 

P*n  P?n 

r?6i 

P3n 

PBn  Pbn‘ 

By  expanding  Equation  a 

set  of  linear  equations 

is  obtained  which 

mav  be  written  as 

A  P  = 

=  B 

(27) 

-a. 

where  A,  P, 

,  and  B 

are  given  by 

(i-k-j 

>2-i 

2^(1-^) 

T^l-kj) 

Tf 

Tf 

Tf/4 

-(l-k1)k2 

-(ki+k2Tf) 

Tf ( l~k1 ) - (T^k^/2 ) 

Tf 

3Tj:/2 

Tj/2 

A  A, 

-(1-k1 

)k3 

-Tfk3 

p(l-k1)-(Tfk3/2)-l 

0 

pTf 

pT?/2 

k2 

-2k? 

-2Tfk2 

0 

2Tf 

4 

k2k3 

“k3 

~ ( P^2+^fk3 ^ 

0 

p-1 

pTf 

k^ 

K3 

0 

-2k3p 

0 

0 

p?-l 

(28a) 

where  k, ,  k„,  and  k.,  are  the  ist,  Pnd,  and  3rd  components  of  the  vec- 
tor  Kp 
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'221 

-kl°E 

In 

in 

-kikfl 

in 

(28b)  and  B  = 

-klk3°E 

in 

,2  2 
~k2°£ 

in 

"k2k3cr^ 

.in. 

,2  2 
rk3°£: 

(28c) 


2  2 

Once  p,  Tf,  ou,  and  ov  (i.e.,  the  model  parameters)  are  speci¬ 
fied,  the  steady-state  gain  vector  and  hence  the  matrix  A  are  both 
uniquely  determined.  Bv  specifying  Or,  the  vector  B  is  also  deter¬ 
mined,  so  that  one  may  easily  solve  for  P  using  a  computer. 


A  computer  program  (named  VPCALC)  was  written  to  solve  for  P 
for  a  wide  variety  of  model  parameters.  A  listing  of  VPCALC  is  given 
in  Appendix  I.  The  main  computation  in  the  program  is  accomplished 
in  subroutine  VPGAUS,  which  forms  the  matrix  A  and  solves  for  P  by 
Gaussian  el iminatioiji  followed  by  back  substitution  8,9,10  .  A  sol¬ 
ution  by  finding  A-*  and  using  the  equation 

P  =  A  1  B 

was  not  persued  because  Gaussian  elimination  followed  by  back  sub¬ 
stitution  is  considerably  more  economical  in  terms  of  computer  time 
than  matrix  inversion. 

Figure  7  (solid  lines)  shows  values  of  /p7 /or ,  the  normalized 
jitter,  obtained  from  VPGAUS  as  a  function  of  panel  R*, 

where  R*  =  ay\2/o  ■  (29) 

For  comparison, j  P^/o  as  given  by  Kalman  theory  for  the  process  and 
observation  models  isvplotted  (dotted  lines)  versus  the  same  param¬ 
eters.  The  latter  results  were  obtained  by  iterating  Equation  (5) 
to  convergence  .  The  standard  deviation  of  the  timing  error  due  to 
measurement  noise  alone  is  observed  to  be  always  less  than  the  standard 
deviation  of  the  timing  error  due  to  the  modeled  process.  Also,  as 
R*  and  p  increase  the  two  curves  tend  to  converge.  Large  R*  and  p 
result  from  models  where  accelerations  are  highly  correlated  and 


iterations  were  halted  when  all  elements  of  P(k)  agreed  to  nine  sig¬ 
nificant  figures  with  the  corresponding  elements  in  P(k-l). 
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measurement  noise  (a  )  dominates  the  maneuver  noise  (a..).  Under  these 

circumstances,  measurement  noise  is  the  predominant  source  of  tracking 
jitter.  Thus,  the  error  standard  deviation  due  to  noise  alone  ap¬ 
proaches  that  due  to  the  model. 

It  should  be  noted  that  the  normalizations  used  in  Figure  7  (i.e., 
R*  and  p)  were  checked  repeatedly  by  running  modified  versions  of 
VPCALC.  „Many  different  versions  were  run  using  different  values  for 
Tp  X,  cr,  and  cr.  In  each  case,  the  resulting  data  agrees  with  that 

presented  here  fJr  corresponding  values  of  R*  and  p. 


B.  Noise  Response  by  Computer  Simulation 


One  method  for  checking  the  validity  of  the  preceding  derivation 
is  Monte-Carlo  simulation.  Here,  the  algorithm  given  by  Equation 
(3)  is  implemented  directly  on  a  computer  with  the  observation  x(k) 
supplied  by  a  random  number  (noise)  generator.  The  state  vector  z(k) 
is  constrained  to  zero  for  all  the  iterations  (i.e.,  all  k).  Since 
the  error  covariance  and  the  predictor  covariance  are  equal  under 
this  constraint,  the  tracking  error  variance  (P,  )  may  be  approximated 
by  the  sample  variance  of  the  1st  component  of  tne  prediction.  That 
is,  for  sufficiently  large  N 


N  0 

l  l?]n{k)|2 

k=l  Jn 
N 


(30) 


where  N  is  the  number  of  iterations.  N  must  be  chosen  large  enough 
to  make  thp  initial  transient  effects  negligible. 

The  computer  program  VSIM,  listed  in  Appendix  II,  performs  the 
above  procedure  for  a  wide  variety  of  model  parameters.  As  the  mod¬ 
eled  acceleration  becomes  more  correlated  (i.e.,  as  p  increases)  the 
steadv-state  time  constant  of  the  corresponding  filter  increases  so 
N  must  be  increased.  Table  1  shows  the  number  of  iterations  used 
in  VSIM  as  a  function  of  p.  The  data  from  the  simulation  is  plotted 
with  circles  on  Figure  7  and  is  virtually  identical  with  that  obtained 
in  Section  A. 

C.  Z-Transform  Approach 

While  Equations  (27)  and  (28)  provide  a  convenient  means  for 
numerically  computing  the  prediction  error  variance,  they  do  not  pro¬ 
vide  a  simple  analytical  formula  where  the  effects  of  each  parameter 
are  observable.  The  desire  for  a  formula  of  the  form 
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N  l 


Table  1 

The  number  of  iterations  per  data  point  in  the  program  VSIM 


0 

Iterations 

.99 

20000 

.999 

40000 

.999 

60000 

.9995 

80000 

.999 

i  oooon 

P-j  =  fl k1 ,k2,k^,Tf,p) 


(311 


led  to  vet  a  third  approach,  Z-transform  analysis.  From  Equation 
(3)  and  the  model  (Equations  (17)1,  a  signal  flow  graph  of  the  Kalman 


Dredictor  may  be  drawn, 
using  Mason's  ru 1 eC 1 53 , 


This  flowgraph  is  shown  in  Figure  8.  By 
the  transfer  function  between  the  input  and 


j  is  found  to  be: 


k,  +  A  Z  1  +  B  Z  2 

H(Z )  =  - - r - a - 

1  +  C  rl  +  0  l~c  +  E  1~ 


02) 


where 


A  =  -k^(l+p)  + 

k2^f  +  k3Tf/'2 

(33a) 

B  =  kjp  -kpT^p 

*  k3TP2 

(33b) 

C  =  k-j  -  p  -  2 

(33c) 

D  =  A  +  2p+l 

(33d) 

E  =  B-p 

(33e) 
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Figure  8.  Signal  flow  graph  of  the  Kalman  Predictor. 
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From  Z-transform  theorytll],  if  x ( k )  has  an  autocorrelation  sequence 
$  (k)  with  transform  $  (z),  the  transformed  autocorrelation  sequence 
oPthe  output  (zj)  is  gPen  by 

4*  *  (Z)  =  H(Z)H(Z'1)$XX(Z)  (34) 

For  the  case  at  hand,  we  have  an  input  of  white  noise,  so  that 

*xx<Z>  •  4  <35> 

and 

♦zjZjfZ)  =  H(Z)H(Z-1)cr|  (36) 

Thus 


k2+A2+B2+(AkjMB)(Z+Z"VkjB(Z2+Z‘2) 
l+C2+D2+E2+(C+CD+DE)fZ+Z'1 WD+CE)  (Z2+Z"2)+E(Z3+Z*3) 

(37) 

Note  that  4>~~(Z)  is  the  Z-transform  of  a  symmetric  autocorrelation 
function,  s6zit  is  a  two-sided  transform.  We  would  like  to  employ 
the  initial  value  theorerrtll]  to  determine  *  (k)|,  which  is  the 

variance  of  the  output.  However,  this  theorim  only  applies  to  a  one 
sided  transform  F(Z)  for  which  f(k)  is  zero  for  negative  k.  In  order 
to  apply  this  theorem,  we  will  split  4>~~(Z)  into  two  parts,  one  cor¬ 
responding  to  non-negative k  and  the  otfier  to  non-positive  k.  That 
is. 


*7-;m 


-j1  (Z)  =  6+(Z)  +  G~(Z)  (38) 

+ 

where  G  (Z)  is  the  transform  of  fk)  for  k  0 

% 

and  G  (Z)  is  the  transform  of  -jr  (k)  for  k  <:  0 

Note  that  at  k=0,  G*(Z)  and  G’(Z)  eacfi  contain  one  half  the  total 
contribution  of  4>~~. 

From  the  Drooerties  of  Z-transformLl?],  if  f(k)  has  the  Z-trans¬ 
form  c(z),  then  ff-kl  has  the  Z  transform  F(Z_1).  Since  —  ( k )  is 
an  even  function, 


(cb^fk)) 

(<b~~(k)). 
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(39) 
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=  {2*(V+CW+DX+EY)+(CV+0W+EX+W+CW+DY)(Z+Z'1)  (47) 

°E 

+  (DV+EW+X+CY) (Z2+Z'2)+(EV+Y) (Z3+Z-3) } 

t  {(1+CZ_1+DZ'2+EZ"3)  (1+CZ+D7.2+EZ3)} 

By  equating  coefficients  of  Z  in  the  numerators  of  Equations  (37) 
and  M 7),  one  obtains 


EV+Y  =  0 

( 48a) 

DV+EW+X+CY  =  k,B 

(48b) 

CV+DW+EX+W+CS+DY  =  Akj+AB 

(48c) 

2(V+CW+DX+EY)  =  kx  +  A2  +  B2 

(48d) 

Note  that  one  could  include  coefficients  for  both  lower  and  higher 
powers  of  Z  in  the  expression  for  Num(Z)  (Equation  (44)).  However, 
after  cross  multiplying  and  equating  coefficients,  one  would  find  that 
the  coefficients  of  these  other  powers  of  Z  are  zero.  Application 
of  the  initial  value  theorem  to  the  expresison  for  G+(Z)  (37)  yields: 


4-do 


=  v 


k=0 


Thus,  we  will  solve  the  set  of  Equations  (48)  for  V.  From  (48a), 


Y  =  -EV 

(49) 

and  by  substitution  for  Y  the  remaining  equations 

can  be  written 

(D-CE)V  +  EW  +  X  =  kjB 

(50a) 

(C-DE)V  +  (0+l)W+(C+E)X  =  Akj+AB 

(50b) 

2(1 -E2) V  +  2CW+2DX  =  k2  +  A2  +  B2 

(50c) 

These  may  be  solved  via  Cramer's  Rule  to  yield 
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^8 

E 

1 

Ak-j+AB 

D+l 

C+E 

k2+A2+B2 

2C 

2D 

D-CE 

E 

1 

C-DE 

D+l 

C+E 

2(1.-E2) 

2C 

2D 

where  |*|  indicates  determinant.  Thus, 

V  =  {(kjB)  (D+l) (20)  +  E(C+E)(k^+A2+B2)  +  2C(Ak1+AB)  (52) 
-(D+l)(k2+A2+B2)-2E(Ak1+AB)D-2k1B  (C+E)C }  + 

{(D-CE) (0+1) (2D)  +  E(C+E)2(1-E2)  +  2(C-DE)C 
-  ( D+l )2(l-E2)-2E(C-DE)D-2( D-CE ) (C+E)C) 


4^  (0) 

Finally,  — ^ —  is  given  by 


<t>~~(0)  *i~(o)  <t>~~(0) 

—  =  ^  ^ 

ct;  at  at 

l  K  £, 


2<J>~~(0) 

T” 


2V 


(53) 


To  check  the  validity  of  Equations  (51)  and  (52)  a  computer  pro¬ 
gram  named  ZSIM  was  written  to  evaluate  them  directly.  A  listing 
of  ZSIM  appears  in  Appendix  III.  For  the  cases  of  p=0.99  and  p=0.995 


the  values  for  J  pin/°r  obtained  from  running  ZSIM  agree  with  those 
presented  in  Figuren7  to  five  significant  figures.  For  the  larger 
values  of  p  (0.999,  0.9995,  and  0.9999)  numerical  instabilities  (i.e., 
over  and/or  under  flows)  were  observed  for  the  larger  values  of  R*. 

The  cause  of  the  instabilities  was  traced  to  the  determinants  in  Eq¬ 
uation  (51).  For  large  R*  and  large  p,  both  of  the  determinants  in 
Equation  (51)  approach  zero.  In  evaluating  the  determinants,  (via 
Equation  (52))  subtraction  of  nearly  equal  terms  to  yield  nearly  zero 
determinants  result  in  a  loss  of  significant  bits  in  the  computation. 
If  the  machine  precision  (i.e.,  word  length)  is  not  sufficient,  all 
of  the  significant  bits  may  be  lost.  The  initial  runs  of  ZSIM  were 
made  in  double  precision  on  a  Datacraft  6024  computer  with  24  bits 
per  word,  resulting  in  approximately  13  decimal  digits  of  precision. 
Many  different  versions  of  ZSIM  were  run,  each  using  different  combi¬ 
nations  of  model  parameters  (o  ,  a;,  T,,  X)  to  generate  R*  and  p. 

None  of  these  were  any  more  orvlesi  successful  than  the  version  given 
in  Appendix  III. 
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As  a  final  check  on  the  cause  of  the  numerical  instabilities 
of  ZSIM,  a  modified  version  was  run  on  an  Amdahl  470  computer  (32 
bit  words)  usinq  extended  (quadruple)  precision.  This  yields  128 
bits  per  variable  or  about  34  decimal  digits  of  precision.  By  com¬ 
puting  the  gains  to  30  significant  figures,  it  was  possible  to  com¬ 
pute  /FT”  /or  to  agree  with  the  data  in  Figure  7  to  at  least 
P  digits;  Thus,  Equation  (52)  appears  to  be  correct  but  must  be  re¬ 
garded  as  numerically  unstable. 
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CHAPTER  IV 

NOISELESS  RESPONSE  OF  THE  KALMAN  PREDICTOR  TO  AIRCRAFT  MANEUVERS 


Recall  that  we  are  using  superposition  to  find  the  Kalman  clock- 
loop  predictor's  response  to  an  input  composed  of  noise  plus  a  deter¬ 
ministic  signal.  Chapter  III  documents  the  Kalman  predictor's  per¬ 
formance  with  white  noise  inputs.  In  this  chapter,  we  solve  the  2nd 
part  of  the  problem  -  finding  the  response  to  deterministic  inputs. 

To  accomplish  this,  it  was  decided  to  simulate  (by  computer) 
actual  maneuvers  rather  than  step  changes  in  delay,  velocity,  or  ac¬ 
celeration.  Although  step  changes  are  much  easier  to  simulate  on 
the  computer,  with  the  possible  exception  of  steps  in  acceleration 
they  are  not  physically  realistic  of  aircraft  maneuvers.  The  maneu¬ 
vers  which  were  simulated  are  shown  in  Figure  9a  and  Figure  10a.  Both 
are  180°  turns.  In  Figure  9a,  the  initial  and  final  flight  paths 
are  along  the  axis  to  the  satellite.  This  is  called  the  axial  180° 

turn.  In  Figure  10a,  the  initial  and  final  flight  paths  are  perpendi¬ 

cular  to  the  axis  to  the  satellite.  This  is  called  the  orthogonal 
180°  turn.  In  both  cases,  the  radius  and  angular  velocity  are  given 
respectively  by  ( see  Reference  14,  Section  8.4) 

R  =  V?/A  f 54) 

and 

P  =  A/V  (55) 

where  A  =  magnitude  of  the  constant,  radially  inward  acceleration 

and  V  =  magnitude  of  the  initial  velocity. 

Of  course,  sj,nce  we  are  simulating  delay  functions,  A  is  given  in  units 
of  chips/sec  and  V  is  given  in  units  of  chips/sec.  Letting  t=0  at 
9=0,  6  is  then  given  by 


0  »  et 


(58) 


For  the  axial  turn,  the  state  vector  z  (i.e.,  (T  ,  f  ,  T  )^) 
is  given  by  (assuming  that  the  satellite  is  at  infinity)  0  0 


Zj(t)  =  T0ft)  =  J 


T  (0)-Yt/c; 

Vn)  jr  sin  yt; 

To(9)  +  Vt/c; 


t  <  0 

n  +  TtV 
°<tlT 
ttV  , 

n t 


(57a) 
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+  RANGE 

y 

—  -c 

0  =  A/ 

— 

DIRECTION 

SATELLITE 


OF  FLIGHT  (a) 


z.{ t)  =  T0  (0)  -  %  sin  0t  ; 

0<t<  »Va 


0  <  t  <  ^v/A 


s  R®/c  sin  5 1  ; 

0  <  t  < 


Figure  9.  The  axial  180"  turn;  (a)flight  path,  (b,c  and  d) 
components  of  the  state  vector. 
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c  T0  (0) 


- if 


/\ 

SATELLITE 


R  =  v%  e  -  A/v 


(t)  =  T0(0)  —  R/c  cosflt  ; 

o<\<*yi 


(t)  =  Rs/c  sin  ; 

0  <  t  < 


=  r®  2/c  cos  0t  ; 

0<  t  <  7r% 


Figure  10.  The  orthogonal  180°  turn;  (a)f light  path,  (b.c  and  d) 
components  of  the  state  vector. 
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2?( tl  =  t) 


-V/C 


t  <  0 


0  <  t 


ttV 
-  7T 


ttV 

T 


<  t 


f  57b  ^ 


0 

,  ra2 

Z3 ( t; )  =  TQ( t)  =  <j  — —  sin  9t 


t  _<  0 

0  <  t  <  ttV/A 


0 


(57a) 


where  T  is  the  delay  from  the  satellite  to  the  coordinate  origin 
of  the  maneuver.  These  are  sketched  in  Figures  9b,  9c,  and  9d,  re¬ 
spectively. 

Similarly,  the  state  vector  z  for  the  orthogonal  turn  is  given 
by 


z  -,f t >  =  TQ(t) 


'V»>  - 1  ; 

t  <  0 

T0m)  - 1 cos  V  1  : 

0  <  t  < 

li0W)  *  §  ; 

M.  <  t 

A  -  1 

jr  (58a) 


z?(t) 


0 


T  (t)  =<  ^  sin  £  t 

0  I  C  V 


0 


t  _<  0 

0  <  t  <IV 

'  A 

"V  . 

7-  <  t 


(58b) 
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z3^  -  v*> =  i 


o  ; 

•2 

R0  A  t 

-  cos  -  t; 

c  v 

0  ; 


t  <  0 


0 -l -T 

TfV  0 

7T  <  9 


(58c) 


These  are  sketched  in  Figures  10b,  10c,  and  lOd,  respectively. 


In  order  to  study  the  behavior  of  the  Kalman  predictor  when  sub¬ 
ject  to  the  delay  inputs  given  by  Equations  (57a)  or  (58a),  a  computer 
program  named  FLYBY  was  written.  After  setting  the  model  parameters 
Tp  p  and  R*  (to  uniquely  determine  the  filter)  and  the  aircraft's 
initial  velocity  and  radial  acceleration  (to  uniquely  determine  the 
turn),  the  program  iterated  the  Kalman  algorithm  (Equation  (3))  with 
a  noiseless  delay  input  using  steady-state  Kalman  ( i . e . ,  Wiener)  gains. 
The  absolute  value  of  the  tracking  error  was  observed  at  each  itera¬ 
tion,  and  after  completion  of  the  maneuver,  the  maximum  absolute  error 
was  written  to  an  output  array.  The  program  was  designed  such  that 
incremental  chanoes  in  model  parameters  and  aircraft  parameters  (in¬ 
itial  velocity  and  radial  acceleration)  could  be  made  in  suitable 
do-loops.  Thus,  the  response  to  a  wide  ranqe  of  inputs  for  a  wide 
variety  of  model  parameters  was  obtained. 


The  structure  of  FLYBY  will  now  be  discussed  in  more  detail. 

A  listing  of  this  program  appears  in  Appendix  IV.  In  order  to  gen¬ 
erate  the  maneuvers,  two  subroutines  were  written.  Subroutines  TURNl 
(axial  turn)  or  TURN2  (orthogonal  turn)  generate  the  state  vector 
z  by  implementing  Equations  (57a, b,  and  c)  or  Equations  (58a, b,  and 
c)  respectively.  Prior  to  beginning  the  maneuver,  the  filter  is  as¬ 
sumed  to  be  tracking  the  delay  function  perfectly.  Subroutines  INIT1 
and  INIT2  initialize  the  predictor  vector  z  (ZWIG)  to  the  proper  values 
corresponding  to  0=0  for  the  axial  and  orthogonal  turns,  respectively. 
The  absolute  tracking  error  at  each  step  (ERR)  along  with  the  maximum 
absolute  tracking  error  thorugh  tne  maneuver  (ERMAX)  is  calculated 
in  subroutine  ERR0R1.  The  predictior  is  performed  in  subroutine  PRED 
using  steady  state  gains.  In  addition  to  generating  the  state  vec¬ 
tor,  the  TURN  subroutines  generate  a  variable  (ICT)  which  indicates 
the  number  of  iterations  elapsed  since  the  end  of  the  maneuver.  An 
IF  statement  terminates  the  simulation  when  ICT  is  greater  than  1000. 
After  the  simulation  ceases,  the  maximum  absolute  tracking  error  for 
the  maneuver  (normalized  to  the  radius  of  the  turn)  is  written  to  the 
array  X0UT.  A  flow  chart  of  this  sequence  is  given  in  Figure  11. 

One  additional  subroutine,  GAINS,  is  used  to  load  the  gain  vector  (GP) 
from  an  arrav  of  gain  vectors  (XGAIN)  corresponding  to  various  model 
parameters. 
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Figure  11.  Flow  chart  for  the  simulation  portion  of  the  computer 
nronram  FI  YBY. 
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The  data  obtained  from  these  programs  is  plotted  in  Figures  12a, 
h,c,d,  and  e  for  the  axial  turn  and  Figures  13a,b,c,d,  and  e  for  the 
orthogonal  turn.  Each  of  these  figures  shows  maximum  tracking  error 
(normalized  to  the  turn  radius)  as  a  function  of  R*  (o  X^/Oy)  for 
a  specific  value  of  p.  Plotted  are  families  of  curvesvfor  various 
angular  velocities.  Consider, .for  example,  an  aircraft  flying  at 
500  mi/hr  making  a  5  g.  (1  g.  =  9.8  m/sec^)  180°  turn:  What  is  the 
tracking-error  for  a  predictor  with  model  parameters  T-  =  .01  sec, 

R*  =  10  ,  and  o  =  .909?  At  40  mchips/sec,  a  1  meter  change  in  range 

generates  a  delay  of  .133A,  so  in  terms  of  delay  500  mi./jjr.  becomes 
V/cA  =  ?o.8  chips/sec  and  5g  becomes  A/cA  =  6.53  chps/sec  .  Thus, 
the  radius  of  the  turn  is  given  by  R/cA  =  136  chips  and  the  angular 
velocity  is  0  =  .00219  rad/sec.  Turning  to  Figure  12c,  corresponding 
to  an  axial  turn  for  6T^  =  .0025  and  R*  =  10  ,  the  normalized  error 

is  4.5  x  10  .  Multiplying  by  the  radius  of  the  turn  (i.e.,  unnor¬ 

malizing)  yields  a  maximum  tracking  error  of  0.06A.  Similarly  for 
the  orthogonal  turn,  the  maximum  tracking  error  is  found  to  be  0.31A. 


In  general,  the  tracking  error  for  the  orthogonal  turn  is  higher 
that  the  error  for  the  axial  turn.  This  is  due  to  the  step  discon¬ 
tinuity  in  acceleration  for  the  orthogonal  turn  (Figure  lOd).  The 
axial  turn  has  no  step  discontinuities.  Also,  the  tracking  error 
for  both  turns  monotonical ly  increases  as  R*  increases.  This  behavior 
is  opposite  to  the  noise  behavior  of  the  filter  presented  in  Figure 

7  where  both  fK  /o  andjp.  /op  monotonical  ly  decrease  as  R*  increases. 
As  R*  becomes  larger  measurement  noise  tends  to  become  the  predominant 
factor.  The  filter  algorithm  compensates  for  this  by  narrowing  the 
filter  bandwidth.  This  results  in  improved  noise  performance  but 
degraded  tracking  of  maneuvers  because  the  narrow  bandwidth  causes 
the  filter  response  to  lag  behind  the  input. 

The  curves  plotted  in  Figures  13c,  d,  and  e  exhibit  behavior 
which  may  appear  strange  at  first  glance.  The  behavior  was  investi¬ 
gated  and  found  to  be  caused  by  the  response  time  of  the  predictor. 

The  responses  of  the  predictor  to  the  orthogonal  turn  is  composed 
of  an  initial  lag  followed  bv  an  overshoot  and  less  significant  ring¬ 
ing.  For  small  R*,  the  response  is  relatively  fast  so  that  the  maxi¬ 
mum  error  occurs  in  initial  lag.  As  R*  becomes  larger  the  filter 
response  slows  and  the  error  in  the  overshoot  becomes  progressively 
more  significant.  The  "bumps"  in  the  curves  of  Figures  13c,  d,  and 
e  results  from  the  occurrance  of  the  maximum  errors  in  the  first 
overshoot.  As  R*  becomes  even  larger,  the  response  of  the  predictor 
becomes  so  slow  that  the  maneuver  is  completed  before  the  initial 
lag  of  the  response  is  completed.  When  this  occurs,  the  errors  in 
the  lag  response  become  progressively  more  significant.  Eventually 
these  errors  again  become  larger  than  those  in  the  overshoot.  The 
regions  of  the  curves  above  the  "bumps"  corresponds  to  the  maximum 
errors  occuring  in  the  initial  lag  response. 
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Figure  12c. 


The  maximum  delay  tracking  error  for  the  axial  turn  as 
a  function  of  model  and  maneuver  parameters ,  P =0.999 


Figure  13c.  The  maximum  delay  tracking  error  for  the  orthogonal  turn 
as  a  function  of  model  and  maneuver  parameters,  0 =0.999 
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R*  =  crv  \Z/cr- 


figure  13e.  The  maximum  delay  tracking  error  for  the  orthogonal  turn 
as  a  function  of  model  and  maneuver  parameters,  cf  0.9999 


CHAPTER  V 

COMBINED  NOISE  AND  MANEUVER  PERFORMANCE 
OF  THE  KALMAN  TRACKING  LOOP 


A.  Introduction 


In  the  previous  two  chapters  the  performance  of  the  Kalman  clock 
loop  in  response  to  noise  inputs  (Chapter  III)  and  delay  inputs  which 
simulate  the  inputs  generated  by  a  maneuvering  terminal  (Chapter  IV) 
was  found.  At  any  given  time  the  total  delay  tracking  error  z,  is 
the  sum  of  the  error  due  to  noise  z,  and  of  the  error  due  to  maneu¬ 
vers  z-,  .  That  is,  l£n 

.1  f 

m 


Z1 


c 


+  Zn 


(59) 


The  noise  component  of  this  sum  is  a  jitter  term  with  standard 

deviation  JP. j  whereas  the  maneuver  term  is  deterministic.  Thus, 

the  error  may  be  considered  jittering  with  standard  deviation /PT 
about  a  bias  term  z^£  .  The  worst  case  behavior  occurs  when  the  Dias 

(maneuver)  term  is  a  Maximum.  The  maximum  delay  tracking  errors  in 
response  to  two  orientations  of  180°  turns  were  found  in  Chapter  IV. 

In  this  chapter,  three  examples  of  Kalman  clock  loop  predictors 
are  presented.  In  each  case,  the  delay  error  standard  deviation  ob¬ 
tained  from  Kalman  theory  (Jp^)  is  compared  with  (1)  the  delay  error 

standard  deviation  due  to  measurement  noise  only  (Jp7  )  and  (2)  the 
maximum  tracking  error  encountered  in  a  5  g.  (1  g  =  3C8  m/sec),  500 
mi/hr.  orthogonal  180°  turn  with  no  measurement  noise.  In  each  ex¬ 
ample  the  measurement  noise  standard  deviation  is  assumed  to  be  .05A. 
This  is  apDroximately  the  same  value  as  the  measurement  noise  in  the 
current  experimental  system. 

B.  Example  1 


In  this  example  model  parameters  identical  to  those  considered 
in  Reference  6,  Example  1,  pp.  111-116  are  used.  These  are  summarized 
in  Table  2.  A  short  description  of  how  these  parameters  were  obtained 
is  given  here.  Values  for  the  code  rate,  the  average  maneuver  rate 
(A),  and  the  measurement  noise  (o  )  were  either  given  or  assumed. 

The  value  of  a'j  was  derived  usingvEquation  (17)  under  the  assumption 
that  Pq  (the  probability  of  zero  acceleration)  is  0.5,  Pm  (the  proba¬ 
bility  of  maximum  acceleration)  is  zero,  and  A  (the  maximum  permitted 


46 


p 

vehicle  acceleration)  is  5  g.  (49  m/sec  ).  The  value  of  was  chosen 

(by  an  iteration  scheme)  to  be  the  largest  value  which  yields/p,  < 
.05A.  It  has  been  shown[7l  that  timinq  jitters  of  greater  than  :05\ 
may  yield  a  significant  degradation  in  bit  error  probability.  The 
complete  set  of  model  parameters  is  given  in  Table  2. 

From  the  set  of  model  parameters  the  derived  parameters  p  and 
R*  are  found  to  be  .9861  and  7.463  10  ,  respectively.  These  values 

may  be  used  with  Figure  7  to  find  the  tracking  error  standard  deviation 

due  to  measurement  noise  only  (7  P ,  ) .  Unfortunately  a  curve  for  XT. 

=  0.014  is  not  plotted.  However,  By  extrapolating  those  values  that 

are  plotted,  it  is  found  that /"F^/oj-  =  .88,  so  that  the  standard 

deviation  of  the  tracking  error  due  to  measurement  noise  only  is  ./F7 
=  . 044A .  ln 


The  values  of  p  and  R*  may  also  be  used  with  Figure  13  to  find 
the  tracking  error  generated  by  a  terminal  performing  an  orthogonal 
180°  turn  at  500  mi/hr.  (232.52  m/sec)  and  5  g.  (49  m/sec^)  radial 
acceleration.  This  turn  has  a  delay  r  li us  of 


^  =  135.95  chips 

(60) 

and  an  angular  velocity  of 

0  =  0.2192  radians/sec. 

(61) 

Thus, 

flTf  =  0.0152  radian. 

(62) 

Using  Figures  13a,  13b,  and  13c  and  extrapolating  to  obtain  a  value 
for  each  figure  for  0T.  =  .0152,  then  extrapolating  these  three  values 
to  obtain  a  value  for  AT^  =  0.014  yields 


Pmax 

Rad 


1.0  *10 


(63) 


Unnormalizing,  the  maximum  tracking  error  obtained  is  approximately 


.136  chip. 


(64) 
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This  tracking  error  may  be  regarded  as  a  bias  in  the  prediction 
due  to  the  maneuver. 

The  initial  selection  of  Tf  for  this  example  was  based  upon  ob¬ 
taining  a  tolerable  clock  loop  jitter  of  .05A.  In  selecting  this 
value  as  tolerable,  an  implicit  assumption  was  made  that  the  clock 
loop  yields  estimates  that  are  unbiased.  This  is  true  for  a  non-ma¬ 
neuvering  terminal  or  a  terminal  making  random  maneuvers.  However, 
with  deterministic  maneuvers  such  as  the  180°  turn  considered  here, 
biases  in  the  estimate  will  occur.  These  may  or  may  not  be  tolerable, 
depending  upon  the  frequency  of  the  maneuver,  whether  communications 
are  expected  during  the  maneuver,  etc.  In  the  present  case,  a  clock 
loop  operating  with  a  jitter  of  .044A  biased  by  .136A  will  result 
in  a  bit  error  rate  of  five  to  six  times  that  of  a  similar  clock  loop 
operating  unbiased  with  a.jitter  of  0.05A  (assuming  one  operates  with 
a  bit  error  rate  near  lO-4 [7,22]). 


To  minimize  the  bias,  it  may  be  necessary  to  reduce  T^  (i.e., 
to  sample  more  often).  Consider  a  system  with  the  same  model  param¬ 
eters  as  the  previous  system  with  the  exception  of  Tf.  Let  Tf  be 
0.0228  (Example  lb).  Since  the  maneuver  rate  (X)  is  unchanged,  p 
(1-XTf)  becomes  0.995.  Similarly,  OTf  now  becomes  0.005^  Using  Fig¬ 
ure  7,  the  Kalman  prediction  error  standard  deviation  (7  P-^ )  and  the 
prediction  error  standard  deviation  due  to  noise  only  (VP-.  )  are  found 
to  he  .03A  and  ,027a,  respect i vel v.  From  Figure  13c  l corresponding 
top=0.995),  the  maximum  tracking  error  for  the  180°  orthoqonal  turn 
is  found  to  be  0.07A.  The  performance  of  the  clock  loop  filter  with 
this  set  of  model  parameters  is  summarized  in  Table  2  for  both  Tf  = 
.0695  (Example  la)  and  Tf  =  .0228  (Example  lb). 

C.  Example  2 


In  this  example  model  parameters  identical  to  those  considered 
in  Reference  6,  Example  2,  pp.  116-119  are  used.  These  are  summarized 
in  Table  3.  Here,  o  and  X  have  the  same  values  as  Example  1.  The 
code2rate  is  doubledvto  80x10°  chips/sec,  which  doubles  oL  to  5.36A 
/sec.  T^  was  reduced  to  .052  sec.  Given  oc,  o  x  and  A  (chip  duration), 
T,  was  chosen[6J  such  that  the  Kalman  predictor  standard  deviation 
of  the  tracking  error  is 

/Fj  -  . 05A.  (65) 


From  the  model  parameters,  the  derived  parameters  R*  and  p  are  qiven 
h  v 


.  2 

a  >  . 

R*  =  — —  =  3.731  •  10 


(66) 
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Table  ?. 

Model  Darameters  and  performance  data  for  Examples  la  and  lb 


Model  Parameters  1  a  1  b 

Code  rate  C/A)  40x10®  cbips/sec  same 

Measurement  noise  standard  deviation  .OF  chip  same 

f  a  /A) 
v 

Averaqe  maneuver  rate  (X)  0.?  maneuvers/sec  same 

Maneuver  noise  standard  deviation  8  cbips/sec  same 

< Oy/Al 

Sample  interval  CT^l  0.069F  sec  0.0??8  sec 


Derived  Parameters 


Acceleration  correlation  coefficient  0.9861 
(p) 

a  X? 

R*  =  7.463  10'4 


.9QF4 


same 


Maneuver  Parameters  f f or  5  g.,  F00  mi/hr.  turn) 

Delay  radius  of  turn  (Rad/cA)  13F.p^  chips  same 

Angular  velocity  (0)  0.?lQ?  rad/sec  same 

Anqle  subtended  in  one  sample  O.OiF?  radian  .00^0  rad 

interval  fBT^l 

Performance 

Kalman  prediction  delaxerro*'  .06  chip  .O'*  chip 

deviation  fj P1 /A) 

Prediction  delav  error  standard  .044  chip  O.077  chip 

deviation  due  to  measure¬ 
ment  noise  only  fJT1n/A) 

Prediction  tracking  error  due  to  F  q.,  0.176  chip  0.07  chip 

F00  mi/hr.  180 

orthogonal  turn  fe  /A) 

J  max 
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and 


p  =  1  -  XTf  =  0.9896. 

The  5  g.  (49  m./sec'),  500  mi/hr.  (223.52  m/sec)  turn  has  a  ra¬ 
dius  (Rad)  of  1019.6  meters.  This  generates  a  change  in  delay  of 

^  =  271.90  chips.  (67) 

The  anqular  velocity  (6)  for  the  turn  is  0.2192  radian/sec  so 
the  anqle  subtended  in  one  sample  interval  is  given  bv 

6Tf  =  .0114  rad.  (68) 

Prom  Fiqure  "7,  the  normalized  standard  deviation  of  the  error 
due  to  measurement  noise  only  (J  P.  /a  1  is  0.88  (for  0=0.09.  R*=3.73 
•in"’).  Thus,  the  unnormalized  jitter  is 


JPln  =  .044A. 

From  Figure  13a  (corresponding  to  =0.991  the  normalized  maximum. track¬ 
ing  error  (e  /Rad)  for  §T  =.0114  is  extrapolated  to  be  6.4*10  . 

The  maximum  kicking  error  is  found  to  be 

p  =  0.174A. 
max 

A  bias  in  the  local  clock  of  this  magnitude  will  result  in  a  degra¬ 
dation  in  the  bit  error  probability  (P^)  by  at  least  afactor  of  ten 
for  systems  operating  unbiased  with  P^smaller  than  10~^[7,22].  As  in 
Example  1,  this  may  or  may  not  be  tolerable.  Table  3  presents  a  short 
summary  of  this  example. 

D.  Example  3 


This  example  illustrates  that  the  data  of  Chapters  two  and  three 
may  be  used  not  only  to  evaluate  the  performance  of  a  Kalman  predictor, 
but  also  to  select  a  set  of  model  parameters  to  achieve  some  desired 
predictor  performance  criteria.  This  may  be  necessary  because  a  de¬ 
signer  may  not  have  enough  aircraft  data  to  select  the  statistical 
parameters  o^  and  \  with  confidence.  This  may  also  be  desirable  in 
that  the  designer  may  want  to  deliberately  skew  the  model  parameters 
(to  reduce  the  errors  in  tracking  of  maneuvers  while  minimizing  pos¬ 
sible  degradations  in  jitter  performance,  for  example).  The  selec¬ 
tion  of  model  parameters  based  on  the  data  of  Chapters  two  and  three 
is  possible  because  the  data  in  these  chapters  is  independent  of  the 
accuracy  of  the  parameters.  That  is,  given  a  set  of  model  parameters, 
the  performance  of  the  resulting  Kalman  predictor  will  be  as  described 
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Table  3 

Model  parameters  and  performance  data 


Model  Parameters 
Code  rate  n/Al 

Measurement  noise  standard  deviation  (a  /A1 
Average  maneuver  rate  (X) 

Manuever  noise  standard  deviation  (cij/ A) 
Sample  interval  (T^l 

Derived  Parameters 

Acceleration  correlation  coefficient  fp) 


Maneuver  Parameters  (for  5  g.,  500  mi/hr.  turn) 

Delay  radius  of  turn  (Rad/cA) 

Angular  velocity  (©) 

Angle  subtended  in  one  sample  interval  (0Tf) 
Performance 


Kalman  prediction  delay  error  standard 
deviation  (fp^/A) 

Prediction  delay  error  standard  deviation 
due  to  measurement  noise  only  f/P^/A) 

Prediction  tracking  error  due  to  5  g. 


500  mi /hr. 

(e  /A) 
max 


1P0  orthoronal  turn 


for  Example  2 

80x10^  chips/sec 
0.05  chip 

0.20  manuevers/sec 
5.36  chips/sec 
0.05?  sec 

.0Q896 
3.73  10"4 


271.89  chips 
0.2192  rad/sec 
0.0114  rad 

0.05  chip 
0.044  chip 
0.174  chip 
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by  the  figures  of  these  chapters  regardless  of  whether  the  model  ac¬ 
curately  reflects  the  physical  situation.  This  is  not  true  of  the 
Kalman  covariance  matrix  ( P ) ,  which  is  accurate  only  if  the  model 
and  model  parameters  accurately  reflect  the  physical  situation. 

In  Example  la,  given  o  ,  Oj,  and  X,  T,  was  chosen  to  be  .0695 
se&.  in  order  to  yield  a  KaYman  prediction'error  standard  deviation 
(J  P, )  of  .05A.  It  was  observed  that  the  maximum  tracking  error  through 
theJ'5  g.,  500  mi/hr.  180°  turn  was  0.136A  and  that  the  jitter  due 
to  measurement  noise  only  was  .044A.  In  Example  lb  Tf  was  reduced 
while  holding  the  other  model  parameters  constant  in  order  to  reduce 
the  tracking  error  due  to  the  maneuver.  This  resulted  in  a  maximum 
tracking  error  of  .07A  and  a  jitter  due  to  measurement  noise  of  0.0P7A. 
Consider  now  the  following  question:  While  holding  the  sample  inter¬ 
val  (Tf)  and  the  measurement  noise  standard  deviation  constant  at 
the  same  values  as  those  considered  in  Example  lb  ( ,0??8  sec  and  . OS  A, 
respectively),  is  it  possible  to  select  o~  and  X  such  that  the  result- 
ino  filter  tracks  the  180°  turn  with  an  error  less  than  .03  and  has 
a  jitter  due  to  measurement  noise  of  less  than  .OAA? 

To  answer  this  ouestion,  a  value  for  p  land  hence  Xl  will  be 
assumed.  Figures  13  Ta,h,c,d,  or  e)  and  7  will  be  examined  to  find 
a  value  for  R*  such  that  both  trackinci  and  noise  criteria  can  be  met. 

Tf  such  a  value  is  found,  then  the  remaining  oarameter  Ta^T  can  be 
easily  found,  since  a  X  and  R*  have  been  specified.  Tf  no  such 
value  for  R*  is  found,  a  different  value  for  p  wi 11  be  assumed  and 
the  process  repeated.  This  procedure  will  continue  until  a  set  of 
model  parameters  is  found  or  until  all  of  the  plotted  values  for 
have  been  examined  and  shown  to  vield  negative  results. 

Let  p  equal  0.99.  From  Eauation  P3)  X  is  0.439  maneuvers/sec. 

From  Fioure  7,  values  for  R*  are  desired  such  that  P,  /a  <  0.8  fi.e., 

_ _  -  In  v- 

Jp,  .04  when  o  =0.5).  This  criterion  can  be  met  if  R*  is  qreater 
than  5.56  10"^Y  From  Figure  L3a  (corresponding  to  P=0.99),  the  maxi¬ 
mum  tracking  error  for  the  orthogonal  turn  (Rad/cA  =  135.95  chips, 
Tf=0.005)  is  less  than  0.03A  (i.e.,  e  /Rad  £  2.2  •  10  .1  if  R*  is 
less  than  1.0  10”'.  Thus,  a  useable  r ange  for  R*  is 

5.56  10"4  <  R*  <  1.0  •  TO'3  .  (69) 

-4  2 

Choosing  R*  to  be  7.5  •  10  results  in  a  value  for  o^.  of  !2.8?A/sec  . 

From  Figure  7,  the  jitter  due  to  measurement  noise  only  is 


while  the  Kalman  prediction  delay  tracking  jitter  is 


JTX  =  .044A.  (711 

From  Figure  1.1a  the  maximum  tracking  error  due  to  the  180°  orthogo¬ 
nal  turn  is  0.037A.  The  model  parameters  and  performance  data  for 
this  example  are  summarized  in  Table  4. 

Thus,  while  sampling  at  the  same  rate  and  with  the  same  energy 
to  noise  density  ratio  in  the  clock  pulse  (i.e.,  the  same  o  )  as  in 
example  lb,  the  predictor  considered  here  will  track  the  deYay  from 
a  180°  turn  with  much  less  error  (.027A)  than  the  predictor  of  example 
lb  ( .07A) .  The  jitter  due  to  measurement  noise  is  somewhat  larger 
here  (.037A)  than  the  previous  example  (.027A),  but  still  very  rea¬ 
sonable. 

Of  course,  in  skewing  the  model  parameters  to  optimize  tracking 
performance,  the  relation  between  the  model  parameters  and  an  actual 
aircraft  may  be  lost.  That  is,  if  the  model  parameters  were  chosen 
based  upon  the  physical  characteristics  of  an  aircraft  (e.g.,  the 
maximum  acceleration  the  aircraft  is  able  to  withstand),  they  may 
no  longer  be  realistic  with  respect  to  those  same  characteristics. 

For  example,  in  the  case  considered  here  the  value  of  12.82A/secz 
for  Oy  corresponds  to  approximately  9.81  g.,  which  is  beyond  the  cap¬ 
abilities  of  most  present  day  aircraft.  This,  however,  does  not  in¬ 
dicate  problems  with  either  Kalman  theory  on  the  chosen  model.  The 
Kalman  algorithms  are  minimum  mean  (i.e.,  average)  sguare  error  al¬ 
gorithms.  This  study  considers  the  performance  with  respect  to  two 
different  performance  criteria,  the  response  due  to  measurement  noise 
(jitter  performance)  and  the  tracking  errors  resulting  from  180°  turns. 
The  model  parameters  which  are  necessary  to  provide  good  performance 
with  respect  to  these  two  criteria  may  be  quite  different  than  those 
reouired  for  the  minimum  square  error  criterion. 

Obviously,  good  trackinq  performance  for  typical  maneuvers  is 
important.  This  example  shows  that  (1)  good  tracking  performance 
for  typical  maneuvers  is  obtainable  simultaneously  with  good  jitter 
performance  and  (2)  that  the  curves  in  Chapters  III  and  IV  can  be 
used  to  obtain  the  model  parameters  necessary  for  the  implementation 
of  a  predictor  with  this  performance. 
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Table  4 

Model  parameters  and  performance  data  for  Example  3 


Model  Parameters 
Code  rate  n/A' 

Measurement  noise  standard  deviation  (av/A) 
Averaqe  maneuver  rate  ^Al 
Manuever  noise  standard  deviation  (cvj/ Al 
Sample  interval  (T^l 

Derived  Parameters 

Acceleration  correlation  coefficient  (p) 


A0*!0*  cbips/sec 
.ns  chips 

.43%  maneuvers/sec 
1?.R?  chips/sec 
0.0??8  sec. 


O.DQ 

7. 5*10“4 


Manuever  Parameters  for  5  g.,  500  mi/hr.  turn 

Delay  radius  (Rad/cA) 

Angular  velocity  (01 

Angle  subtended  in  one  sample  period  f 0T^ ) 
Performance 


1 35.95  chips 
0.?lh?  rad/sec 
0.005  rad 


Kalman  predictor  delay  error  standard  .044A 

deviation  (JP-j) 

Predictor  delav  error  standard  deviation  .m7A 

due  to  measurement  noise  only 

Predictor  tracking  error  (maximum)  due  .037A 

to  5  g.,  500  mi/hr.  180°  orthogonal 

turn  (em,  ) 
max 
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CHAPTER  VI 

TRANSIENT  PERFORMANCE  OF  THE  KALMAN  FILTER 


A.  Introduction 

In  this  chapter,  some  aspects  of  the  Kalman  and  Wiener  clock 
loop  filters  during  acquisition  or  lock  up  are  examined.  Since  the 
Kalman  filter  is  computationally  more  costly  to  implement  than  the 
Wiener  filter,  knowledge  of  when  it  is  advantageous  to  use  one  over 
the  other  is  needed.  The  important  criteria  for  comparison  are  the 
magnitude  of  the  tracking  errors  (overshoots  and/or  lags)  during  lock 
up  and  the  time  required  for  the  error  to  converge  to  reasonably  small 
values.  Two  computer  programs  were  written,  named  TRANS  and  STEADY, 
to  examine  the  behavior  of  the  Kalman  filter  (TRANSl  and  the  Wiener 
filter  (STEADY)  under  lock  up  conditions.  These  programs  and  the 
data  they  generate  are  examined  in  this  chapter.  Before  discussing 
these  programs,  the  initialization  of  the  Kalman  filter  as  described 
in  Reference  6  is  examined. 


Initialization  of  the  Kalman  TDMA  Clock  Loop  Filter 


In  Reference  6,  it  is  proposed  that  the  initial  estimate  (zf-l)) 
be  obtained  from  two  measurements  prior  to  enabling  the  filter.  The 
observation  models  corresponding  to  the  two  measurements  are 


xf-2)  =  Zj(-2)  +  $-2) 
and 


(7?) 


x(-l)  =  z7 (-1)  +  Ef-1) 

A 

and  the  estimate  zf-l)  is 


f  73) 


( z.(-l)  +  F(-l) 

zJf-D-z^) 

rTr’~ 


(74) 

(75) 

(76) 


Note  that  no  attempt  is  made  to  secure  acceleration  information.  The 
initial  prediction  is  obtained  from  the  filtered  estimate  z(-l)  via 
Equation  (6b)  and  is 

z(0/-l)  =  *z(-l)  .  (77) 
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The  effects  of  noise  and  the  neglected  acceleration  component,  on  the 
initial  prediction  z(0/-l)  will  now  be  analyzed.  Here,  an  acceleration 
with  constant  magnitude  of  T  (chips/sec  )  throughout  the  initiali¬ 
zation  period  as  assumed.  Tnus,  with  this  constraint,  the  state  vec¬ 
tor  z(k)  is  given  by 


( zj_(k)  ^ 

(  Zj(k)  ^ 

z(k)  = 

z2(k) 

= 

z2(k) 

[z3M  j 

^fo  J 

,  TQ  constant  .  (78) 


At  k=-l,  the  1st  two  components  of  the  state  vector  z  are  given  by 


Zj(-l)*  Zj(-2)  +  z2(-2)Tf  +  TQT2f/2  (79a) 

and 

z2(-l)  =  z2(-2)  +  T0Tf  .  (79b) 

From  Equation  (75)  and  Equation  (79a  and  b), 

z2(-l)  =  z2(-l)  -  \  r0Tf  +  (80) 

is  obtained. 

At  k=0,  the  1st  two  components  of  the  state  vector  are 

z-j  (0)  =  Zj(-l)  +  z2(-l)Tf  +  T^/2  (81a) 

z2(0)  =  z2f-l)  +  T0Tf  (81b) 

From  Equation  (77),  the  predicted  values  for  z-^(O)  and  z2(0)  are 

Zj(0/-1)  =  z | ( - 1 )  +  z2(-l)Tf  (82a) 

z2(0/-l)  =  z2(-l)  (82b) 


By  substituting  z.(-l)  from  Equation  (74)  and  z?(-l)  from  Equation 
(80)  into  Equations  (82a)  and  (82b)  we  obtain 

z^OM)  =  Zj(-l)  +  z2(-l)Tf  -  \  ToT2  +  2c(-l)-5(-2) 

which  can  be  rewritten  using  z^(0)  from  Equation  (81a)  as 

Zj(0/-1)  =  zj(0)  -  TQT*  +  2e(-1)  -  c(-2)  (83) 
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Thus,  the  prediction  z,(0/-l)  can  be  considered  biased  by  an  accelera¬ 
tion  term  of  -In  andxcorrupted  by  2  noise  samples,  F(-l)  and 
F(-2).  Similarly7 for  the  velocity  component  of  the  state  vector, 
by  substituting  £,,(-11  from  Equation  (80)  into  Equation  (82b)  we  ob¬ 
tain 

z2(0/-l)  =  z2(-l)  -  \  T0Tf  +  ^lL'-|Il?.)_  (84) 

From  Equation  (81b),  z2(-1)  =  Zp(0)-TQT^t  so  Equation  f 84)  becomes 
z?(0/-l)  =  z2(0)  -  |  T0Tf  +  (85) 

Thus,  the  velocity  prediction  is  biased  by  an  acceleration  term  of 
-  3/2  T  Tf  and  corrupted  by  2  noise  samples  weighted  by  1/Tf.  For 
PN  code  rates  of  40 * 1 0°  chips/sec  and  T^=.01  sec  (quite  large,  con¬ 
sidering  the  chip  rate),  the  bias  errors' introduced  into  the  delay 
and  velocity  predictions  are  1.3*10  2  chips  and  2*10~'5  chips'sec, 
respectively,  per  6  (1  6  =  9.8  m/sec  )  of  acceleration.  Both  of  these 
are  indeed  negligible. 

From  Equations  (83)  and  (85),  the  errors  in  the  delay  and  ve¬ 
locity  components  of  the  estimates  due  to  noise  are  respectively 


Z.  (0)  =  25(-l)-C(-2) 
n 

(86a) 

*2«  <°>  ■ 

(86b) 

n  'f 


The  error  variances  associated  with  these  terms  are  (since  the  noise 
samples  are  uncorrelated) 


Ef(z,  (0))?\  =  5ol 
en  ^ 

(87a) 

and  , 

?  2ar 

Ef<V°»  ’  ■  / 

'f 

(87b) 

with  corresponding  standard  deviations  of 

Oj  *  J Ef(ikj0))2}  =  2.24  o^ 

(88a) 
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and 


(88b) 


Equations  (88a)  and  (88b)  are  important  statistics  because  they  give 
information  about  the  magnitude  of  the  errors.  Assuming  that  the 
errors  are  Gaussian  distributed,  the  prediction  (zf0/-1))  will  be  with 
in  one  standard  deviation  of  the  correct  value  68%  of  the  time  and 
within  1.65  standard  deviations  90%  of  the  time[13].  For  example, 
if  O£=0.05A  (a  reasonable  figure  for  the  present  system),  then  the 
errors  in  the  initial  delay  prediction  z,(0)  would  be  less  than  .185A 
90%  of  the  time.  Simi larly,~for  =  0.01  sec.,  the  errors  in  the 
initial  velocity  prediction  z?(0)  would  be  less  than  7.05A/sec  (or 
.07A/Tjr)  90%  of  the  time.  The  programs  TRANS  and  STEADY  simulate 
the  behavior  of  Kalman  and  Wiener  filters  during  lock-up  with  initial 
errors  of  these  magnitudes. 


C.  The  Steady  State  Kalman  Filter  Program 


The  program  STEADY  (Appendix  V)  is  structured  to  make  extensive 
use  of  subroutines.  Because  the  input  signal  to  the  Kalman  filters 
is  the  difference  between  the  time  bases  of  the  local  and  received 
clock  signals,  initialization  with  an  error  in  any  of  the  components 
of  prediction  z (O'- 1 )  is  equivalent  to  a  step  change  of  opposite  sign 
in  the  corresponding  component  of  state  vector  z(0). 

Three  subroutines  were  written  which  generate  the  state  vector 
simulating  step  changes  in  delay  (T  ),  velocity  (t  ),  or  acceleration 
(T  ).  These  are  names  STEP,  RAMP,  and  PARAB,  respectively.  The  names 
originate  because  a  step  change  in  velocity  produces  a- ramp  in  delay 
and  step  change  in  acceleration  produces  a  parabolic  delay. 

After  the  state  vector  is  generated,  it  is  compared  with  the 
prediction  in  subroutine  ERR0R1  or  ERROR?.  These  subroutines  keep 
track  of  the  errors  as  the  filter  converges.  It  was  discovered  that 
for  step  changes  in  delay  (T  ),  the  maximum  error  almost  always  oc¬ 
curs  at  the  first  iteration  and  is  equal  to  the  magnitude  of  the  step. 
When  this  is  not  the  case,  the  maximum  error  occurs  in  the  first  over¬ 
shoot.  Accordingly,  ERROR?  was  written  which  computes  the  magnitude 
of  the  first  overshoot.  This  subroutine  is  used  for  step  changes 
in  delay.  For  steps  in  velocity  and  acceleration,  ERROR!  is  used. 

This  subroutine  computes  the  maximum  error  as  the  filter  converges 
regardless  of  where  it  occurs. 
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After  the  error  is  observed,  two  subroutines  are  used  to  check 
convergence.  The  first,  subroutine  COUNT,  counts  the  number  of  it¬ 
erations  for  the  first  < the  delay)  component  of  the  prediction  to 
converqe  to  within  .OFA  of  the  actual  delay.  This  is  used  to  obtain 
data  concerning  the  convergence  time  of  the  filter.  The  second  sub¬ 
routine.  STOP0,  is  used  to  determine  when  all  three  components  of? 
z(0)  have  converged  to  within  .001  (chips,  chips/sec,  or  chips/sec' 
depending  upon  the  component)  of  the  corresponding  actual  values. 

If  all  three  components  have  converged,  a  variable  ISTOP  is  set  to 
one  which  ceases  the  iterations.  If  not,  subroutine  PRED  is  used 
to  generate  a  new  prediction  by  direct  implementation  of  Equation 
13).  Subroutine  PRED  may  be  used  with  either  time  variable  or  steady- 
state  gains.  It  should  be  noted  that  subroutine  PRED  implements  a 
linear  tracking  characteristic.  The  nonlinear  characteristic  of  the 
correlation  circuitry  could  be  simulated,  but  this  would  prohibit  the 
normalization  of  output  data.  From  this  point,  program  flow  is  de¬ 
pendent  upon  whether  the  filter  has  converged  (by  the  criteria  of 
ST0P2)  or  not.  If  it  has  converged,  then  either  the  number  of  it¬ 
erations  (from  COUNT)  to  converge  or  the  error  (from  ERROR!  or  ERROR?) 
during  the  convergence  process  may  be  written  in  an  output  array. 

If  the  filter  has  not  converged,  then  the  next  state  is  generated 
and  the  process  is  repeated. 

D.  The  Transient  Kalman  Filter  Program 

The  program  TRANS  1 Appendix  VI )  uses  the  same  subroutines  and 
is  organized  in  the  same  manner  as  STEADY.  Three  additional  subrou¬ 
tines  are  also  used  for  functions  unique  to  the  time  variable  (Kalman) 
filter.  _/irst,  it  is  necessary  to  obtain  an  initial  error  covariance 
matrix,  P(0).  In  Ref.  6  it  is  shown  that  this  matrix  should  be 
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Subroutine  SETUP?  implements  Equation  (89)  directly.  Second,  it  is 
necessary  to  calculate  the  predictor  gains.  This  is  accomplished 
in  subroutine  GPCALC  via  direct  implementation  of  Equation  (4).  In 
this  equation,  the  gains  are  calculated  as  a  function  of  the  error 
covariance  matrix,  hence  the  need  for  SETUP2  prior  to  GPCALC.  Third, 
after  the  prediction  is  obtained,  it  is  necessary  to  update  the  error 
covariance  matrix  for  the  next  iteration.  This  is  performed  in  sub¬ 
routine  PWt  by  implementing  Equation  (5). 

E.  Results 


These  computer  proqrams  were  used  to  find  M)  the  errors  in  de¬ 
lay  tracking  z.  during  lock-up  and  (?)  the  number  of  iterations  for 
the  delay  error  to  fall  below  .05  chips  for  a  wide  variety  of  model 
parameters  and  inputs.  The  magnitude  of  the  delay  errors  are  pre¬ 
sented  in  Figures  14,  15,  and  16  for  inputs  of  steos,  ramps,  and  pa¬ 
rabolas,  respectively.  In  each  case,  the  results  are  normalized  to 
the  input  to  make  the  curves  more  general.  That  is,  the  data  for 
the  step  input  is  normalized  to  T  the  ramp  data  to  f  Tf,  and  the 
parabola  data  to  T  T^/2.  In  all  three  figures,  data  forthe  Kalman 
filter  is  presente8  in  solid  lines  and  data  for  the  Weiner  filter 
is  presented  in  dotted  lines. 

While  developing,  and  debugging  of  the  programs,  it  was  deter¬ 
mined  that  the  maximum  errors  in  the  ramp  and  parabola  cases  are  under¬ 
shoots  (or  lags)  in  the  response,  whereas  the  errors  in  the  step  input 
case  are  overshoots.  In  Figure  15  the  data  for  the  Kalman  response 
lies  on  the  R*  axis,  not  below  it. 

The  data  for  convergence  is  presented  in  Figures  17  and  18  as 
a  function  of  model  parameters  for  the  Wiener  and  Kalman  filters. 

Figure  17  shows  the  iterations  to  converge  given  a  step  offset  in 
delav  of  0.70A.  This  corresponds  to  the  90%  confidence  interval  for 
Oj-=0.05A,  Tf=.01.  Figure  18  shows  the  iterations  to  converge  for 
a  step  change  in  delav  rate  of  7  chips/sec.  Aqain,  this  corresponds 
to  the  90%  confidence  interval  for  Or=. 05A  and  T^=.0i  sec.  These 
fiqures  are  specific  for  the  inputs  described  ( tney  are  not  normal¬ 
ized).  However,  the  general  characteristics  may  be  generalized. 

First,  the  most  significant  improvement  of  the  Kalman  filter  over 
the  Wiener  filter  occurs  for  large  R*  and  larqe  p.  This  corresponds 
to  a  relatively  slow  responding  filter  in  the  steady  state.  Second, 
the  number  of  iterations  required  for  converoence  of  the  Kalman  pre¬ 
dictor  seems  to  be  asymptotic  to  13  in  Figure  17  and  31  in  Fiqure 
18.  This  behavior  is  explainable.  In  Reference  6  it  is  shown  that 
the  Kalman  gain  vectors  for  different  model  parameters  all  follow 
the  same  monotonical ly  decreasing  trajectory  until  very  near  this 
final  (steady-state)  value.  Consider  two  filters,  the  first  of  which 
has  a  gain  vector  which  converges  to  (within  1%  of)  steady-state  in 
?00  iterations  and  the  second  of  which  has  a  gain  vector  which  con¬ 
verges  (1%)  in  400  iterations.  Suppose  that  both  filters  are  given 
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the  same  input,  and  the  first  filter  acquires  the  signal  within  100 
iterations.  Under  these  conditions  the  second  filter  will  also  ac¬ 
quire  the  signal  within  100  iterations  because  the  gain  vectors  for 
the  two  filters  are  almost  identical  until  near  the  200th  iteration. 
Furthermore,  any  filter  whose  gains  go  to  steady  state  in  signifi¬ 
cantly  more  than  100  iterations  will  have  similar  convergence  charac¬ 
teristics.  Third,  for  the  slower  filters,  the  difference  between 
the  Kalman  and  Wiener  filters  is  dramatic.  Consider  R*  =  10"  ,  p=.999, 
and  T,=.0J.  Figure  17  shows  that  the  Kalman  filter  requires  13  it¬ 
erations  to  converge  (.13  sec  at  T.=.01  seel  from  an  initial  delay 
offset  of  . 2A  whereas  the  Wiener  filter  requires  150  iterations  (1.5 
seel.  Figure  14,  when  interpreted  for  T  =  0.2A,  shows  that  the  1st 
overshoot  is  .1A  for  the  Kalman  case  and°0.07A  for  the  Wiener  case. 
Thus,  the  filter  remains  in  the  linear  portion  of  the  correlation 
circuitry  throughout  lock-up  (see  Figure  3).  Similarly,  for  a  ve¬ 
locity  offset  of  7A/sec  (Figure  18),  the  Kalman  filter  requires  31 
iterations  to  converge  while  the  Wiener  filter  required  over  700  it¬ 
erations  to  converge.  From  Figure  15  the  maximum  errors  during  con¬ 
vergence  are  .07  (Kalman)  and  1.4A  (Wiener).  Since  the  correlation 
circuitry  is  linear  only  over  a  ±^A  range,  the  Wiener  filter  would 
probably  not  lock  under  these  conditions. 


Magnitude  of  the  maximum  predictor  tracking  error  in 
response  to  step  changes  in  delay  acceleration  of 
magnitude  To. 


ITERATIONS 


R  *  =  crv  X2/crf 


Figure  18.  The  number  of  iterations  for  the  predictor  tracking  error 

to  diminish  to  less  than  0.05Agiven  a  step  change  in  delay 
velocity  of  7^ /sec. 


CHAPTER  VII 
CONCLUSIONS 


The  steady  state  and  transient  performance  of  a  Kalman  filter/pre¬ 
dictor  implementation  of  the  receive  clock  loop  in  a  TDMA  satellite 
system  has  been  analyzed.  The  steady-state  response  of  the  filter 
to  measurement  noise  was  determined  by  state-variable  methods  and 
verified  by  Monte-Carlo  simulation  and  z-transform  techniques.  The 
steady-state  responses  of  the  filter  to  delay  changes  caused  by  a 
terminal  performing  180°  tu^ns  with  two  different  orientations  with 
respect  to  the  satellite  was  determined  by  computer  simulation. 

Examples  are  qiven  which  show  how  the  steadv  state  responses 
may  be  used  in  the  evaluation  and  design  of  the  clock  loop  filter. 

The  examples  given  show  that  the  errors  in  response  to  terminals  per- 
formino  180°  turns  are  somewhat  larqer  than  those  one  would  intui¬ 
tively  expect  from  the  data  qiven  bv  the  Kalman  error  covariance  matrix. 
The  Kalman  algorithms,  of  course,  minimize  the  average  square  error, 
not  the  error  due  to  a  single  deterministic  maneuver.  Nevertheless, 
the  response  to  typical  maneuvers  is  important.  The  examples  show 
how  it  is  possible  to  deliberately  skew  the  model  parameters  of  the 
Kalman  algorithms  to  improve  the  response  to  typical  deterministic 
maneuvers  with  little  degradation  in  noise  performance  (jitter). 

The  transient  performance  of  the  Kalman  clock  loop  filter  and 
the  corresponding  Wiener  clock  loop  filter  was  analyzed  by  computer 
simulation.  The  tracking  errors  caused  by  errors  in  initialization 
for  both  Wiener  and  Kalman  predictors  were  determined  as  a  function 
of  the  model  parameters  of  interest.  The  iterations  to  converge  were 
also  determined  as  a  function  of  model  parameters  for  two  typical 
acquisitions.  The  results  show  that  the  Kalman  filter  has  much  better 
transient  performance  than  the  corresponding  Wiener  filter  when  the 
derived  model  parameters  R*  and  p  are  both  relatively  large. 
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Appendix  I 

The  Computer  Program  VPCALC 


(_****•»*****■**»■**  PkOUPAK  VPCALC  *************** 

OPTIONS  DP 

u  i  1 1 cs i u;  pi; (6 > ,  peel 6 > , x < 3 > .  xhat< 3  > , op ( 3  > , xnoi  om > 

DI  YENS  I  CL  XOUT (4,4  1 ,5) ,V(6) 
lit  AL  LAMUAt  ti) 

CO'u.‘. Ofl/PHI  /rd'O  ,T(4  ) 

CU OK/  5  E'1'2/ S  I GV ,  V  AR  V ,  SI  GW .  V  A R tv .  S I GA  .  V  A P  A 
DATA  LAI.BDA/  I .  . .  5,  .  I  . .  .’5  . .  0 1  / 

SI  GV  =  . 00 

VAIiV=SICV**2 

CALL  DH/SSN 

'1(  I  )=D.l'l 

T ( 2 ) =0.  i.'.JC!  I 

TO)  ='•).  l  C.’OC  I 

iU)=i'.uriT  ciM 

r  ACTOK=  KL  **('.  I 

"0  v,M  ILA'  =  I ,  j 

f>;  !0=  I  -LA  "EDA  ( I  LA  M )  *T  ( I  > 

r(  =  U>  •  *  *— E>/  r  ACT  01) 

F.'O  bl.NJ  IPT=I,4I 
li=u*rACTOi) 

si  c/=s  i  c.-v*i.a.’.:dda  ( i  la;.'>  **p/r- 

V  AKA=S  I  (.  A**2 

VA..1  =VAI>A*(  I  .-RH0**2  ) 

SI  0:  =SOtiT(  VA.ii//) 

I r ( I  FT .  i:0 .  I  )  CALL  SETUP  (Pa1  ) 

C* *•*■•*’••  I TEi; ATE  TO  STEADY  STATE 
lELAOS 

00  l.liJ  1  =  1,1  coco 
CALL  Pl.’l  (PL, POO) 

CA  LL  CHHC,:  I  (  PI. .  PI-.O,  I  FLAG ) 
ir  (  It-LAC  .EO.  I  )G0  TO  I  O 

I  Ci-  COI.jLTUE 

>  *>  vnr:,  *  1: rt  i.  *  x  -k  7  A  ***  x  *  *★  ieie  * 

I I  \j  ,'u  1  2*1  1  =  1  «  j 

x: 'at  n 

141  0l.LT  l  NU'ti 

CALL  ;'PCALC(  f  :.,0P) 

CALL  VG/L'S  (CP,  V  ) 

Xt.’CT  (1,1  i'T  ,  I  LAP  )  =p 

aoiit  (2 ,  i  pi  ,  i la:  )  «sorn ( ( 1  )/v  akv  > 

X TUT (  3 ,  IPT,ILA!'.)«SOin(v(  !  >/VAWV  ) 

AoLT  (4,  ITT.  ILA  :>=■.  .0 

CALL  ASS  1  GL<  6K>ATAS<i ,5'U  t  *3S,2) 

.  11,’  lt:(P).-'T.T 


h  rr-i  cl 'EOT  1  (  p  1 

. ;  i I L  ?..(o),  P..r  (r;) (i > ) 
1C  1=1 ,6 


Ks  V  -  X.  -  -.A  . 


'■  nun  pkaoUga JUiS 

**caj  lOuDG 
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be  IKPUin.En.K.lCO  TO  II 

b'l  P(  I>=(  AlSIPHfl  >-PWO<  I )  ))/APS(Ph(  I  )  ) 

be  JO  Tu  It. 

bv  ii  P(i)=ij. 

cu  ii.  cuiniMit 

O  I  L»G  lb  1=1,6 

o2  lb  Ir(P( I ) .Of. It .**—  I  I . )GO  TO  50 

o_  I ELAG= I 

c-  PbTUnl! 

Ob  bt.  1  FLAG=U 

to  HETUHN 

c7  El  ID 

08  C1-***’  4*  *-*  **■**».•****************  ********** 

cv  SUliKOUTI f.E  DECGMPCNN.A.IIL) 

7 1  DU’E.'JSICN  A(0,6) ,UL(6,6) ,SCALES(6> ,IPS(6) 

7  1  CO!. KOM  .IPS 

7 1  N=I!N 

7  3  C 

74  C  I  Nil  IAI.IZE  I  PS. UL,  ANT  SCALES 

7  b  DO  5  I  =  I  ,  N 

7  o  IPS(I)  =  I 

7  7  HOKJiMKsC-.il 

7  6  DO  2  J=l  ,M 

7  v  UL (I  ,J  )=A<  I , J > 

hW  IK  ( HOl.NI.  K-AI-'S  ( lil.  ( I  .J  :>)  I  .2,2 

t  i  i  ho;vnmi<=afs(ut.(  i,  jj  ) 

L2  2  COIIT IllUb 

t3  IK(K0I/NKH)  3,4,j 

64  3  SCALD S(  I  )  =  I  .  C/H01.JIMB 

bb  GC  TO  b 

fco  ‘  CALL  Slf G<  1  ) 

r 7  SCALES ( I >  =0. C 

bo  b  CONTI: II :t 

LV  C 

V.  c  GAUSS  I  At:  ELIMINATION  WITH  OAHTIAt.  PIV^TI'T 

VI  N.'.i  I  ='l-  I 

vi  DU  IV  K=  I ,  I.1’-'  I 

V3  UIG=O..W 

V4  DO  II  I =K,  N 

Vb  IP=IPS(I) 

vc  SI 2b= Arts  (l)t.(  IP.K  J  )*SCA>  HS(  IP  ) 

V 7  Ir (SIZE-RIG)  I  I  .  I  1 ,  IT 

VS  it)  L I G=SIZL 

V V  I CXP I  V  =  I 

I. iti  II  CONTI  HUE 

It  I  IKCBIG)  13,  12,13 

It2  12  CALL  SILG( 2 ) 

I  if.  3  GO  1 0  I  / 

104  13  IKU)X?I\i-K)  14,15,  U 

lib  i  4  j'=I  PSCK  ) 

I I  o  I  p S { ti )  =  I  PS  C I IVX P I  V  ) 

107  I  PS(  I  i)XF  I  V  )  =  J 

106  lb  rCPalPS(tC) 

10V  PIVOT«t;L(KP.K) 

II. )  KP  I  =.C*  I 
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‘  .  I  .  G  I  9  I  =  i\P  1  ,  * 

I  1C  IF  =  |P'fc(l) 

I  .  -  =- 1 1 1_ (  I  f-..\)/PI  ,  ■/, 

ii4  i:l  U'v: 

i  1  N  i  C  I  9  C =*\P  1  , N 

i  ‘0  '.I. ( :  )+H,’*uL<:<r.J ) 

o',  v.  i.’-.r.c  lc(  p.  lift-:  .‘‘Ac hi mf  lahoijaoe  coni ng  it  cdm^h/--;- 

Ilf  v.  .’Lt:C  :.:!•]  ^rtCL'lIvCE  rrrICIENT  CODE 

i .  ■»  u.  con  1 1  nut 

I.  .  l'<  CvN'i  I  NN'E 

I  <r  I  ,\>-  =  I ) 

UC  Ir(  jUril  ,!l);  IV.Irt,  lv 

U--  lo  CALL  Lit  GC2) 

1  c. ,  iv  »t  I'Ciit! 


V  *  it  *•*■*★"*  *’.  k'k  -t 


ii'iv  sc.v>:c:.,,;,ui..h,x) 

I  SUti  l!L(o,6),i-(c),X(o),  1PM6) 


I  o  ^  C 

I r=I PS ( I ) 

I--  X(  I  )=.!(  IP) 

I-:  CO  2  1=2,:: 

to.  i  p = ;  r  s  c : ) 

io‘,  I  11  =  1-1 

lot  LU  / 

1  -  l,  o  O  I  o  =  I  ,  I  ’  M 

1‘k  i  SL'“  =  f)L".'+';L  ( I  P,  J )  *X  ( J  ) 

Ul  i  X(  I  )=E(  IP) -SL 

i  AS  C 

i  *  -  i -=  i  ps  c > 

:<(t:)=xa:)/uuip,:i) 

I  No  I;0  A  I  tAC,<  =2  ,  N 

i-t  i=;tpt-iiACK 

in',  c  i  goes  (N-n . 2,i 

I  Nil  I  P=I  PS  (  I  ) 

1  ‘  V  I  P  I  =  I  +  I 

I  PC  hU‘.'  =  ,'.rt 

it i  ih;  j  j =]  p i  ,t: 

IL2  -  SIJ.'f=S'J‘.!+lIL(IP,J)*X(J) 

lb  2  <-  X(  I  )  =  (X<  I  )-SI!‘l )/!!!.( IP,  I  ) 

ISn  RETURN 

i‘-o  END 

I  L  C  C *****  ******************************* 

i  b"i  SOBkOOII t'E  SIt!G(  Il'IHY  ) 

Ibri  p  1  FORMAT!  LNHO’.'ATRI  X  1MTH  ZERO  ROD  IM  DECOMPOSE. 

) 

IbV  12  rOh.>.:Al'(5<HB£i;iGI)LAB  MATRIX  IN  PECONPOSF.  ZERO  DIVIDE  IN  SOLVE 
) 

I  to  12  FGKMATCtNHt  NO  CONVERGENCE  IM  IMPRUV.  "ATRIX  IS  NEAR!*  SP'Gip|./R 
) 

I  (■■  I  NCUT =9 

let  o  LOUT-STAKDARO  OHTPMT  UNIT 

I  to  Go  10  (1, 2, 2), I  I'HY 

ICN  1  ;..iITF(t.Un  ,  11  ) 

i  (.  t>  iJU  ,0  H 


■XlVtS  ^ 


-T 

EU 


practical3 
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i  CO 

i.rtl  i  I: (  NLUT  ,  1  2  ) 

l  G*i 

GO  3  0  Ik. 

1  06 

..!<[  IctKUIT,  13) 

1  ov  lo 

WtTUifl 

1  ~i  L. 

tf.'b 

1  t  1  U*  *irn'A  •kitit'k  **■*>  ■**■*■*  •**■*'*  *★■**★★**■*•*  *•★ 

I  '(2  SUBKCUTI  Ni:  VPCALC(GPtV) 

172  niMfc'l.'SICr:  AtC.o)  ,!!L(ft,ft>  ,0P<3>  ,V(6) ,  R(  6) 

1 1 *•  cOiV'.MOii/pri /«Hu, i  c 4 ) 

I  7  b  CO/. KON/iRl 2/ n  GV  .  V ARV , S 1 0.'.' .  V  ARi'-' ,  SI  GA  ,  V  AR  A 

17c  At  I  ,  I  )  =  (  1  .  -GPt  I  )  )**2-  I  . 

I"/'/  At  I  ,2)=2.*<  I  .-GPt  I  ))*T<  1  ) 

17 d  At  I ,  .">)  =  (  I  .-GPt  I  )  )*Tt2) 

I'/S.  At  I  , 4)  =T ( 2  ) 

At  I  ,  b )  =7  1 3  ) 

Ibl  At  I ,o)=7t4 )/4 . 

I s2  A(2, I )=-( I ,-CPt I  ))*GP(2) 

I c 2  At 2, 2  )=-G:> (  I  )-GPt2  )*T(  I  ) 

lo-  A ( 2 , 3 ) = ( i .-Of  t  1 ) )*Tt  I ) -GPt  2 ) *T  f  2  )/2. 

I  c_  A  ( 2 , 4  )  =7  t  I  ) 

Ice  A(2,b)=;.*T(2)/2. 

Id 7  A (2 , 6)=  7  f 3  )/2 . 

it.c;  A  (  3 ,  I  )=-(  I  .-CPt  I  ) )  *GPt  3 ) 

I  c  V  A  ( 3 , 2  !  =-C  P  1 3  )  *7  (  I  ) 

IV.  A  (  3 , 3  )  =  <  I  .-GPt  I  )  )*RHO-GPf3)*T(2)/2  I  . 

I  V  I  A  ( 3 , 4  )  =(, . 

I  V2  At  3,  i> )  =i.KO*7  t  I  ) 

IS-;  A  (  3 ,  c ) =h  1 :0*T  1 2  )  /2  . 

IV*  A ( 4 , I )=CF ( 2 )**2 

I  V 3  A ( 4 , 2 ) =-2 . *GP ( 2 ) 

IVc  A ( 4 ,3 ) =-2 . *GP( 2 ) *7 ( I ) 

IV  7  A(4,4)  =  !:. 

1  Vo  At  4, 5) =2 .*T(  I ) 

! VV  A ( 4 , 6 ) =T ( 2 > 

2t  L  A  <  5i ,  I  )=GP(  2  )*GPt  3) 

2t.  I  At  5, 2  )  =-GP  (3  ) 

21.2  A  (  0 , 3 )  =-l-.HC*CP  f  2  )-T  (  I  )*GP(3) 

21.3  A  ( 5 , 4 )  =i. . 

A  t  b,  s>)  =1  PC- 1  . 

2;:b  A  (  3 , 6  )  =I.HO*T  1 1  ) 

2t  P  A  (  ft,  I  )=CP  t  3  )**2 

207  At  c,  2  )  =t. . 

2i.c  A  (  6, 3)=-2  .  *RHO*GP(  3 ) 

2c  V  A(6,4)=t. . 

u  a  t  o ,  j )  =t  . 

211  A ( 6 , 6 ) =hHO**2- I . 

2  12  Pt  I  ;=-VAl.V*GP(  I  )  **2 

2  13  Bt  2  )  =-VA!;V*GP(  I  )*GP(  2  ) 

2  14  iit3)=-V/kV*GPt  I  )*GP(3) 

2  1b  B  t  4 ) =-VA  RV*OP ( 2 ) **2 

2  1 1  bt  5 )  =-VA  !;vt*'GPt  2  )  *G?(  3  ) 

217  l;  1 6 )  =- VA 1.  V*GP  t  3 )  -**2 

21d  /i,'i=6 

2  IV  ..'ALL  UECC;\t>(LU'A,UD 

22 it  CA LL  LOLVt;  (  LI  .  UL ,  b ,  V  ) 


3 


£ 
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r 


1 


2ts. 
23  I 


2  3<* 


>  I 


c 

2  +  / 


4  +  <: 
2  ‘  *y 


RET'JhN 

2ND 

SUbRl.'UVI  f'H  VCAOS  ((>:», OP) 

1 ’  i ns  ter  pv  (o),GP(3) 

CU'-'.U.N/PHI  /KH0.K4  ) 

eO.,.!.,OK/iET2/SIGV  , V ARV, SI  GO , V  ARE , SI GA , V  AK A 
P=  I  ./(  P't  (  I  )+V ARV  ) 

GP(  1  )  =  (  Pi  t  I  )+Pi.(2)*T<  t  )  +  PI< ( 3  )*T ( 2  )  /2  . )  *P 

UP (2  )=(F>.(2)*P!.(3)*V<  I  ))*p 

GP(3)=Pt.(3)*Ki!0*P 

hi:  TURN 

t'M 

***»  *-.;*x,;?.***-Jk*******.**'********* 

SUniCIJTIf.'E  ?!•  1  (PC/,  PNC) 

>_1  31' l  ROUTINE  CALCULATES  THE  CNR  STEP  PREDICTOR  ERROR 

Utnnno/AKIAIiCLT.  IT  USES  THE  PRESENT  CNF  ST  FP  PREDICT  nr  V  AP  I  a  ".JFS 
«.****» aNII  l-.bTL’hUS  MTU  UPDATED  ONE  STEP  VARI  Af.CFSC  P.\ )  A’")  THE 
o  il;  STEP  VARIANCES (PM)) 
o r ‘ E  :sic::  piuoi.pawo) 

CU:-.'Cii/rMI /iil'j,  i  (4  ) 

Cl, I  .51  :/i  Hi 2/ £  I UV  ,  V  Ai< V ,  SI  0  v .  V ARE . S [ GA  , ,  At:  A 
CA  f.'.V  A=  l  . 

Du  I  [=1,6 
p..o(  i  )=;  i  ( i  > 
it  CONTINUE 

p«(p.:c;  1 1  i+VAiiv ) 

A=E..(5(  I  )+T(  I  )*PaG(2)+T(2  )*P.'/0(  :  )  /?  . 

e=ga:  ,va»p;-.o(2)+T( i  )*pko( :■) 

P.<(  I  )=P>  C(  I  )+2  ,*T(  I  )*P«0(2  )+T<  2)*(  P>«:<  3)+p.'/CU  )  )+T(3  >*?aT.(5) 

3  +7(4)  *f'i'.G(  6 )  /4 .  -  A**2/P 

p.. (2  ;=r./  f  ;.'.a*f.:o<2  )+T <  t  >*p;to<  :)+GA"><A*r  <  i  )*p;,pi  4)+(  i .  ♦oavu/?  . ) 
i.  »’T( 2  ) *:N'  C( E  )  +  T (  3) *PI  C(6  )/2. -A*R/P 
»  3  )=L:  n*(Pi  0(3  )  +  T  (  1  )  (5  )+T(2 )  *P'.‘.'o  (rt  )/2 .  -A  +  pn-H  3  )/0) 

•;  (4  )=GAD  M A**2*P ISC  (4  )+2 .  *GA.'“‘ A*T  (  I  )*o.-'0<  5  >  *7  <?  >*P  “  <  o ) -P-+*5/t: 

:  ..  (s  )=!.T  C*  ((>/  A*UCO  ( b  )+T(  I  )  *p;.C(6)-p+  M  0(  3  )/’  > 

Pm  ( o  )  =  ( Pr  0(o  0(  3)**2/P)*:f!'0**2+V  AT 

RET  DM.' 


*  *  *.;•*  **  -a  »•+*-*?■  **  *+*  *******  •**■**  + 

SUbt-DUTH":  SETUP  (PA  5 
.3 : e:.s i c  n  Pi-(o) 

COT  V  )/rHI/,<I  C,T(4) 

CO.  J:  Vi  E.  2/5 IGV  . At;  V  ,  Si  G.:  .VAR  l-  , S I GA , v  AR  A 
'  ..(  I  )  =  ’.'  /  r.  / 

(3  )  =  . ;  ;/T  (  i  ) 

’..CO  )  = 

!  ,  (4  )  =  vr|.A*'i  (2  )/4.  +2  ,*V ARV/T ( 2  ) 

1  1.(0)=.:;  C  +  V/;.A*T(  I  1/2. 

:  •  .Do)*,/,:  / 


t„  a'  t  A  *  v  7  •;  ^  ^  *  *  >  -A- 5V  ■*  •>  ★ 


EBAJWCIJiUS 

'  •  ’  J.  V>.  ••  ‘  •  —  ,  ,; >£ 


f*- 


.  >  0  R+/Vrf 


I 
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Appendix  II 

The  Computer  Program  VSIM 


*  C 

7 


o 

~7 

3  J 

At; 

4  1 
2 


A  C 

4;^ 
4  V 


.VAHL.SIOA.VARA 

’1/ 


2  options  m> 

include  ‘Pssr.sysv 

-  DMKISK?:  XC)  .XHATC2)  .GAIN’S  (7 .4  1 ,5)  ,X0UT<4,4I  ,5) 

5  REAL  LA.’-.F:DAC  L* ) 

c  OOiMiV.Ofl/PHI  /HI  0.1(4  ) 

COL  ;..0N/ t  E » 2/ 5 1 UV , V AH  V , SI  OH 
d  DATA  LAI  rNA/l  . ,  .5, .  1  ,.M5.. 

v  T<D=.m 

1 1,  1  (2)=T(  I  )**2 

1  i  1 1 3  )  =’j  (2  )*Y(  I  ) 

12  !( 4 ) =T<2 ) **2 

I-  CALL  LE/SSI! 

I-  CALL  ASLIGC(5Hfi!-Pl.  I.5H4I('.3S,  I) 

id  RE  AIM I ) CAINS 

If  C  L  OS  li  I 

t r'ACYOh=  I  (' .  **i  .  I 
I.;  CALL  XORI'HKX) 

IV  5=1. 

2  2  A.'i=<). 

2  I  I X=7227 I 

22  1)0  9 .4.1  ila;i=i,5 

c:  r:u..:=ii.A.-*2('iVi) 

24  I<ii0=  I  .C-LALEDAM  I  LA.V)*T(  I  ) 

2d  N  = 1 0 • **-5/E ACTOR 

IX)  1  PT=  I  ,4  1 

k=i.*E  AC'i  (.  m‘ 

C*****  iOt:  7;iE  SI  Mill.  ATI  ON 
CALL  XChpOl  (XIIAT) 

F.W-.'.i' 

S  in  A  1;  =  «; .  (‘ 

c  ALL  A£l  I(5t:(  ^•'Sr:Er,',r)K4l:l5S.  I  ) 
l-.nAiM  I  )  I  X 
CLOSE  I 
no  -j.;o  i  =  i 

CALL  ItKEEr  (  X  .  Al'AT.SliM ,  S'.«  EAU  ) 

CALL  CAL:S(  IX,  S,  AN. X NO  I  ) 

X I  =(,' . 25*  XNOI 

CALI.  (OAL.St  I ,  IP  t,  I  La.:  ),XHAT,X,XI) 

up  conti  :;ue 

CALL  Ail  I<;N(4:iSFEI),5E4  l.:;s.  I  ) 

UilTEt  l)IX 
CLOSE  1 
SU.*'=n  Li  *  /filJ  A- 
S...EAN=SI.=AIVi  HU 
XCIITf  I  ,  I  PT  ,  1 1. A  "• )  =R 

XCL i  (2,1  £>(  ,II.A  ')=i;A[‘  S(7,IPT,II.AM) 

XCU'i  C  ,  I  £*1  .  ILAL)=?ORT(St|"  I/.05 
xol  ;  (4,  i  ,  I i.a::)  =si'Ea:i 

CALL  ASLI0L(6I")ATAA5,5H4  l!l.-S,2  ) 

..«•  i .  .a?  ) ;  * 'A  it 

CLOSE  7 
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i 


<  X?  ir  r  *  x  ■*  ■x  it  ★  ★  ★ 


i :  ;•;  \ i  t :-. ) 

i,l'  :•  -.it.:;  xc  ) 

.\  (  I  )  =.  . 
a'2  )  =.: . 

a(  j)=.;. 

i.'fc  i  UhN 


k  *x*5.  *A  *  ’Ar  x  MicJrt*  *  *  *  *  *  x  H  x  -A-*****-*** 


I’LHn^uu::^  :</::'  :;ux,  tY.YrL) 
!  Y=1  a*  I  cc<- v 
!  r  l 1  Y )  -i ,  c ,  o 
o  I  Y=  I.  +  C c  I 

■:  Vt- r  =  i  v 

'■:  r  L  =  YrL>  .  i  !  V.Wj  ;--0 


c/usim >..  ?,  \  i,v  > 


L<.  V.'  I  «  I  ,  -Sf» 

call  .< a; ::::.!( u,  iv.ri 

u=i  Y 
A  =  A  +  V 

V  =  ( A  /A  .  -  I L’  .  !*S+A  M 

hLTLAIi 

ti!:i  • 


Sl.T.rtO’JYli®  P ii LD (  GP ,  X !JAT ,  X ,  X I  ) 
i Li :s icy  op c j ) . xkat c 3 > , x < ? ) 

cu;.v:u;/n;i/iii-'o.T(4) 

xii:cv«:u  t  )->.i:at(  i  >+xi 

CALL  SUf  Lull <  C  ) 

XiiAL  ( 1  )  =  ;  :iA'i  ( 1  )+t<  1 ) *xhat< ?  i  +t  <2  >/?  .*>:pat<  ?)  +  ;=<  1  )*:•■  p'ov 
XiiA’i  (2  )  =  >.i  A  L(P  >  +  T<  I  )  *.<r;ATC-  )  +:p(  ? >  *X  I  L'CV 
LHAI  C-.  )=!-.liC*>,MA".  (3  )  +  r:,(.3  )*>  I  vnv 
CALL  SiIPL:,;<(-3) 

•AbUitc: 


*  *-.V  •?: \V  **  -kn  **  **  x^Hif***^*!  **J.* 

Lui’hCUi  I .’II  E-KEEP  CX.XI'AT.SIIV,  5MEAV) 
Or'LLLIL::  XC)  ,XHAT (3) 

Li(i,=  ::t  I ) “ x I ■  a *1  (  i ) 

S,  =  :A;:+L!-.i) 


f--k  ★★  ***★  *Vr  X*-  ★  ★ 


f  ;v*w  * 


Appendix  III 

The  Computer  Program  ZSIM 


..•  *  x*x  xx  1  i;  - }  '*f\  A.‘  .  Zl  I  w  ***tst*xr****.-* 

LPVIuNS 

I : .vJL '  II : "  cSof  ,  3  r  0  + 

; ■  i : 'c  1 1  r  o ( 7 , i . 3  ),>(■-,  ^  i . 3 ) ,  c  p  ( 3 ) 
r-.AL  la: 

c  .  /  !/.!■.,  > 

:  A : a  la;.  '■  -'t/  i ... 3,  .  i , . .’-J ,  .PI / 

it  i  )  =  ..-i 
i  <  a  >  =  i  ( ; )  **;’ 

; j )  - 1  ( 1 5 •*!(?) 

.  (-  >=i(j)  ’■*?. 

CALI.  /:>  I •  >: ' (  LHO' I  ,  5H4  I  P3f  .  1  ) 

•.  i:  A!  ’  (  I  )  C 
CL-.-..-  I 

r  AC  .  1  C  •  +*  ' .  I 

v,  n.;"=i  ,-j 

n  H  .=  I  .-LA  /•(  I  LA"  5*7  (  I  > 

. .  =  I  .  *+— .  /r  AC7 1  t 

call  oi: i 
call  (3) 

CALL  LCf  P-w.-t-'.) 

ro  ;;>:=i.-:i 

K=.*;*r  AC7CI; 

op c  i  )•>•::(  i ,  IP7.  it. a;:) 

,  )=oi.\  nv. ,  m.ac  > 

•JP(3)*G<3,  IP7,  I  LA  I ) 

CALL  ZLL3  ACS') 

a<  i ,  ii>;, ila" )=;; 

X<2,  I  PC,  ILA."  )=0(  7.  ITT,  I  LAI ) 

L  (3,  IPT.ILA’  5  =  SOIi7  (ACS) 
a(4,  ir-r,  ila; 

CO.'.1 1  I.T.Ii: 

C<  ■ .  i  1 1  .UP. 

CALL  AOnc;  (ol-C'ATACI  ,51141. :35,I  ) 

Li.  iv:.<  i  :x 
CLOSI:  I 
CALL  bin 


C  l  ire:  Z Fi:2 (OP  , A LS ) 

-CLU;  GP(3) 

-  ;:i 1 3  ,)ii'.',3.r.nc 

o,v(4) 

;|’(  I  )*(  ii!  !C  +  I  .  )  +  GP(2  )*T(  I  )  +OP(  3 )  *T  (7  5/2. 

■(!)-■  c«)0-rr <2  )*7(1 ) *k I '0+0° ( 3)*T<2>/2. 

'(  I  0 

i- :  . +*  i’o 
-’■IiC 

=  (  I  .  +  ro*2  .*n+p*(C+F)*2.*<  I .  -H**2  ) 

mo-;  *•-•  )*;-.*c-i  .*(  1  .+r))*2.*(  1  .-f**.’  ) 

(  0-0*  r. )  *2  .  *L-  ( L-C*i:  )*  ( C+E  )  *2  .  *C 

=  !!*():  (  I  )*  (  I  .  +  0*2  .*!'+  F*  (C+F  )*(GP(  1)  **2  +  A**2+F**2  > 

*  (  A  *(  L  (  I  )  +  A*P  1*2  .*C-  I  .*<  1  .  +  11)  *(  OP  (  I  )**2  +  A**2  +  r;**2  I 
i  A*nr  ( 1  )+/+h)  *2 . +r;-2*op  ( 1  )*cc+E  )*2.  *c 


l 
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Appendix  IV 

The  Computer  Program  FLYBY 


V 
I . 

I  i 

I  il 

I  3 
i «. 
t  t 
l  o 


2  c 


3  i 


(_************7**  p [!GOi< a  i  t- 1  ypy  **************** 

cru-j.-js  up 

DI.'-lf.'L-K’;  2(5),Z:-.'IOC3).OP(3)  ,X0UT(7,4I  .5) 
1. 1  •Ir.'JijR.’'  XG/IM3.4I  ,5) 

KfcAL  LA,'‘R.')A(5) 

(_•  jv.  j.;/n:i/i<i-'o,i(4) 

cu;  ..•w;:/J-2T3//OC.  VKL, THED.HAn.PI 

!.'A  1  A  la;.  (  DA/  .5,.4.2  5,i,.H5,lJ.  025  .0.005/ 

GALL  ASS  I  ON  ( WI’-'ESSG,  5H4  I  l'.3S,6) 

GALL  DLASSf.' 

PI =3 . M \5v2 A 5355/o 

CALL  A  S3 I  UK(  LHOC'OOS  ,  5114  I  0JS,  I) 

h'LAIJ  f  I  )  /.CUT 

CLOSE  l 

i)0  1/  1  =  1,3 

r.o  i  ,:  j  =  i ,  4  i 

r/j  10  l  = i ,5 

i.  agahi i , j. .•;)=;•  out <  i  ,j.k) 

vo  2/  i= : ,  7 

n.j  2. '  j=i.4| 

I’C  .<*1,5 

«...  >..;i 'id. j. <>=i  . 

■;  ( i  )=..i2 
. 1 2  J  = i C  J  )  *  *2 
T  C  3 )  =T  (  i  )**.-. 

.  \  '■ )  -  \  (  I  I 

r  AGi  .'!•■=  i  P.  **v  .  I 

:;o  v  .  ;  !_a "  =  i ,  5 

«:'C=  I  .-LA  l>;  /(H.  ) *T(  I  ) 

I ::  I.=  3 .; . 

ACC=P . * i 3 . 

:;o  y,v  i :  i  i , 7 
I  ;;:<P=ILi::;:-! 

ACC=Aa.V2 . 

3:!L:i=ACC/.  I  ! 

a  AO—  /.;L**2/’  X 

.,=  10.  **-  5/f-ACTOK 

'.  w  1T=I  I 

call  a r;  i  ,]'( 4f=HrA;,5P4  i:i3S.  I ) 

i-ki  n.ir.’H,  ipv,  it. a  i 


...  =  a*r  AC  ,  ( 

CALL  I • : I T ?  < X V IO.RiLAX.ICT) 

call  c.Ai  ::so.c.Aiii(  i .  i  pt.  I  la;.  )  .op) 

::c  ;  =  i ,  i < 

i  ,;ci  ) 

If  I  C  .  . 1 . .  ! .  PC  )  ( 10  75c 

C.M.I  -■'!  •'•.)(,•(  1)  .XV  !!'•(  1)  ,HH;  .-WAX ) 

lall  c.i  ,  or  > 


ALL  /  ■>:.  i  i  c  :-LVc-  ,  ll’4  i  .  i ) 

I  .  ; :  i  ;  >  v.  i  ;  ; 


::v 
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be  CLOSE  I 

b7  CALL  LX  I T 

bo  END 

V  C*  ***■*■■** -a  *-*•■***  y.  ***■/**■*  ***★  ***** 
ot  v,»***>S,JBhOirn::i:  ERROR  I  COMPUTES  THE  WAX  I  MU,"  ERROR  (OVERSHOOT) 
cl  SUBROUTINE  ERRORI  (Z.ZUIG.ERR.ERMAX) 

o2  c  nK~ Z~ /.  L I  ( j 

ob  ABLR'K=AL  S(  LKK ) 

L*.  If  (E;  .1AX.LT.  ABLER)  ERMAX=AFERR 

cb  RETURN 

co  emu 

C  #  L’**;iA1''******1',rA**  AJli******^^**** 

CO  C* **** SUFiktiU I'l  NK  PLED  COMPUTES  TI  E  OLE  STEP  KALMAN  PREDICTION 
OV  l>****USU!0  ECU  AH  CM  3.  HERE  THE  VARIABLE  X  IS  THE 
7 C  (.»***»  lilPU I  AM)  lit  < : V  13  THE  I Ni.PVAI I CH. 

7  1  SURRUUl'l  r>:  PI.IHH  v.z:  10, OP) 

72  -  I  Mo "L I  <  L’  Zl. IG( 3  ) ,  CP  (3  ) 

7c  REAL  Ii.’NCV 

7 cci.:.o;!/! hi /j;i  o,T(c) 

7b  i;,lio-/=;;-z:.'I(;(  i  ) 

7c  z  ’io ( i  io ( i  )+/.;.■  io (2  )*t<  i  )+zito(3>*T(2 )/?. +op(  i  )*innov 

7  7  /(  10(2)  =2  1  10(2  )+Zi'  IG  (3  >*j(  I  )  +0P(  2  )*I  "NOV 

7  ri  ZM(iCr)«i.1KI«;  .'I  <)  (  >  )+(:P  ( 3  )  *  I  HI.'OV 

7 V  RETURN 

M. 

>r  I  C*  ****  A*  *i.  rtS**.  *****  **  *  **  rVr’r  ***** 

E2  C»*>«*%SUb...!o,T.VE  CAINS  TEA'  SELLS  T*L  CAINS  L'RO'I  THE  ARRAY 


t: 

XC-A.r  Tl  THE  V ECTOR 

CP 

LA 

LohnOUilllr:  GAIUSOIC 

AIL,  ,'rp  ) 

C  5 

UI MEHSICN  XG/IC(3), 

G  f-  <  V 

Ld 

op(  i  )=;n.  a  i  n  ( i ) 

tf/ 

GP  ( 2  )  =  XL  A  I  !■'(  2  ) 

to 

GP(3)=XGA!L'(.-) 

bv 

RLTUi::: 

s  <. 

tin) 

Sr  1 

C,'*  *  *  ■*•*,  **  *  *  w  -A  *  v  *  >;  -A  *  *  *.v.<  -.Hr*..- 

*y.?c*** 

SUUnOUTI  HE  7  DEM  (/., 

J. I  XXX) 

V- 

.OI.'iENSUN  ZC) 

S-4 

uu;.;v.::/pi:i/ei  o.rm 

Vb 

CUV.,"' ’.'./S.ET3/ACC,  VET 

.IHF-'.RAD.'M 

Vo 

I  r-  <  i  > a;<  .(;:: .  i  k;<:  to 

20 

ELI  7' IT. 


I  I 2 


I  1  ‘* 

J  i  c 
i  1 

I  i  c 
i  I  V 
Id. 
12  I 


2  o 
24 


j  I 

Ju 


4C 

4 ", 


^4 

CO 

< 

tiV 
C  it/ 

o  I 

c  j 
6- 
oc 


Z(  1  )=«:(*.  >**  1  *  ;. 

Kciw.-!,'. 

L  ‘A  x  x  >  *  ’.i  x  x  «  n  ■*  >  *  S  ^  •*  •*•  *  *  *  *  *:V  * 

(,****»  xl!l^  ar.CMU.  CC^xiiih  j>!K  X  VECICrf  |:U..i;";  Pit 
t-LV'Ji/.  PrrdOl)  Ar’PTT  Tiir  ."-AI'E UV Fl( . 

*;  t.  I  A  ,\  -- 1  A  /  /  ■*"  I 

z ( i )  =z ( i  )+::(z)*r(  i  ; 

r » b  i  u  .■ : : : 

‘4*  ****  **  **•„*  r.  w  **  t*x*sKS*S**-^* 

SOKKt  ItilTI  (ZriG.KHUAX.ICT) 

ui.M=:;iKt!  z;  iotj  j 
cc : :/i  ••: t'3/agc ,  vh i..  v it';.* ,  ■•.  ad  , p  i 
z.-u;c  i  )=u. 

Z.\IU(2)=-VI:L 

Z  i‘>  I  O'  (  4  )  = ’  • 

!:K..'AX=.'. 

ICi  =  : 
ziiT!'.;; 

4  7f  >  *  *★',***  ■*  V  *  V*  *  >'  f  -V  A*  X  4  7T  * 

o oi: .>4 l/n  r  U  1  ?'I  12  (2 1'  I  < Z.<  '.<  ,  IC7 ) 

iT'HNSitr.  zi. !';<;•  ) 
c  j.  • ( ):i/f ■  r  i  'j/t  c;;,  / 2  L . T  n  t; ) ,  ■<  t  o , «  i 

i-  ■>  I  0  (  I  )=-K.Mi 
.2  i  >  I  G  ( 2  )  -  i  . 

Z/.IGC  j)-f. 
fcTHAXsP . 

ICT=  ■ 

KfcT  ’.li<  -I 
h.'iD 

[J*  %‘Ar* v  *  “A  *Ar SI  **x  •>* *•/*  Vix*** 

SULMOUVINS  TUi;?'2  (Z,J  ,  I  XXX  ) 

zc-) 

cc:'..;o:;/r;  \/\\\ 

co;  S =  i  3//cc,  vEL,'Hi?.:),;<Ar:.r>t 
IF  (I  XI.X.OE  . !  TO  HO¬ 
TKEY  I  )  *(  J-l  ) 

I  r  (7: \r.7l  .GT.  PI  >00  TO  I  .1 
U  I  )=-(.'  /I  *001  (VP  FT  A) 

Z(2 )  =+ka;.*ttf:.*:,I'-  ni'ETA) 
yjji  =+,:/ [;*: *:t 2**2* ccoc  <  theta  j  > 

RETUiTi 

t>*x*>*X**iH  ★  »  Si  -A  VfVr-A  ieici.-  -k inkle ’Jr  kick  v- 

C**x**TMIS  GECT!  ON*  Cu'«r*UTE$  TUH  >  VECTOR  P'l.-r^  ?"!• 
c*****Vi<a:-si  nor:  neon  .‘lAftcvi-.-ar  a  aio't  to  stra ictr; 
l>:  I  xxx=*i 

ALpHA=  :>  ST  A— f- 1 
ii  .ib=A!.r!iA/Ti!.;;j 
Z(  j)=J.( 

2(2  )=ti. 

Z(  I  )=.<AL 
kI-Tuio; 

(J>  *  *  *  A  r.  *>★  *A  ■*  ■A*  ie-J;  'hi;  ■***•>•> 

(,«**** Tf-IS  StCTICIi  CO." PITIES  THE  X  VECTOR  ri'MIIHG  THE 


•  t  --TL I« 


riuxciiCAW* 


yp  ■* 


r'UCHT. 
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to  c»****rLYOl.'i  r  ri,<  1 01  Arlini  THK  '  :i-Mvhk  . 

C  i  i,  l.  I  X  A  A  =  I  A  X  X  +  I 

Oii  7(1)  i 

cv  hclUh!; 

'<  i.  ii.x'i: 

I  t  lj***^n*X***V  ?.> 

7z  ?u[;:<ou7i;:i£  ii:i.::3<z,j,ic7> 

i  ~  i  1 1 1  i: I  f  b  I  L ! !  Z  (  ~  ) 

/<:  cor:.<a;/i  ::i/k’Ki),7(4  ) 

i-j  ;j;i/tH73//-cc, v'::L,7;i=n,nAi},i>i 

7  c  1 1-  (I  C  i'.IT.  I  )  GO  70  1  : 

7  7  (  1  3  -V  ( J  —  1  ) 

7  8  IM'lfifcY/  .(37.  P!  )  HO  ,r  |„i 

7  v  Z<  j)=+U/i'*i!IEU*-!r2--r(C('S(7H;:"i  A)-l  .  ) 

to  Z(2)=Z(2)+Z(j)*T( I ) 

ti  z<  i  )=z<  i  )+Z(;>)*7(2)+:/(;)+;o>)/2. 

t:Z  IHtJ.cO.  I)Z(I  )=-.iA!l 

t3  ir(J.c"(:.  I  )Z( 2 )=■•:. 

t4-  Ht'.'UiO: 

re  it.  IOIO'i+l 

t  c  Z  < I >  =Z (  I ) 

t  7  2  <  2  )  =Z  ( X  ; 

tt;  2  (  3  )  =0 . 

tv  khTlOo; 

v  7  !-;;!]  i 
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Appendix  V 

The  Computer  Program  STEADY 


1  c*  **** **********  PROGRAM  STEADY  *************** 

2  GPTIO.'iS  DP 

-  DIMENSION  2C)  .ZiVlGC-),GD(3)  XOUTI 4, 4  I  ,  5  ) ,  XGA I  »H  7, 4  I  ,5  > 

4  HEAL  LAI' ED  A  <  5) 

b  CCM>:  ON/PHI  /UKO.TC  4 ) 

c  DATA  LAf.  BD A/ 1 . . . 5 , . I , . 05  , .  D I  / 

'/  CALL  IJEASSN 

O  CALL  ASSIGN!  5IIGOMO  I  ,SH4  11*35,  I ) 

v  kEALM  I  >\GAIN 

I 2  CLOSE  I 

11  rm  =  .oi 

12  1 ( 2)=T( I )**2 

IS  i  (  3  )  =T  (  I  )  **3 

14  T(4)*i( I )**4 

ib  i-AC i  Ct\=  I !).  **(’,  I 

It  DC  y  EC  I  CNT:<I.=  1 , 3 

I 7  !  nj  y  •  >;  1 1.A  ‘  =  1,5 

I  o  KHO*  I  .-LA  *r.n/<  I  LAM  )*T(  I  i 

1  *  r  =  !  • '■ .  **— *j/h  A  C  *  OP 

s,  do  ipy=i,4i 

2  i  r!  =  R*r  ACT  CD 

1.  .ALL  GA]"S(XOAI'!(  I  .IPT.ILAM)  ,GP) 

2. -.  CALL  ii:T.»(Z,2r'HO.cK.V/X) 

24  D;)  7;';:  1  =  1, bt'iNW 

Lb  JO  '.':)  (  It', 2W  ,3.1 )  ,  I  CN7RL 

X  il  CALL  STTPtZ) 

27  CALL  :-hi.OH2(Z(  I  )  ,  Z T.  I  0(  I  )  .cPR.EH.MAX) 

2  0  CO  TO  be 

2y  *i.  CAl.L  i.A..  ? (  2  ,7(  I  )  ) 

CO  TO  4., 

i  I  t.  CALI.  PAI- A  ri  ( Z  ,T  (  I  !) 

22  CALI.  EIT-.CE  I  <ZM  )  ,Zi'.'lG(  I  )  .E'RR.ERMAX) 

...  to  CALL  STT  -X  Hl<R ,  HR  MAX,  I  STOP ) 
o  CALL  Pi;ED(Z(  1  )  .Z.'IIG.GP) 

2o  (.-( I  STOP  .SO.  I  )GO  TO  l()i! 

4  0  7  i  -1 1  CONVINCE 

-7  Tv.  XOL  1  (  I  ,  i  °Y , 1  LA". )  =l< 

TO  (o  ',  70, HO)  ,  ICNIRI. 

-  y  i.i.  X'JOT  (2, 1  RT,  I  I.Af.')  =E'R.MAX 
- -.  G(;  TO  y( 

4  i  7..  ;<UuY  ( ,  I  :-7  ,  !  LAM )  =ER.MAX/T (  I  ) 

47  DC  7  7  op 

‘  -  c._  ,JJ'  !T  (  A  ,  I  "T  ,  I L  AM  )  =cif.MAX/7  (2  ) 

“•  v.  CALL  A SS  I G:: ( 6USK POiJT . 5N4  I D3S ,  I  ) 

..  n  I  i  2  (  I  i  X  OUT 

-7  > .  .■  <:„P.  I  Hit 

4,-.  y:  1  c Nit 

‘V  VS.  Co'.  71.'. Ul- 

CM.L  EXIT 

•.  i  No 

•j  ,'  ★>;:**.  ■>  *A  **  :h'»  *>  *  *  -Ar*  *★*"*•  ★ 

...  o>  ..a  *»  NONE  Co  i  I  OH  STEP  GENERATES  A  (UNIT)  SV'>. 

.  ..  S  U:  i-  ACT  I  Nr:  EIKP(Z) 

‘  :  i  Or  J!l 


■  ^ 


KLl** 


3- 
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r 


bo  Z  (  I  )  =  l  . 

b7  Z(2)=U. 

b8  Z  (  3 )  =t) . 

bv  RETURN 

CO  END 

O  i  C*  ***^»  *★★■*'  *  ***  ★★*★★★  ->*★ 

62  c*****f UBKOUTJ t'E  .'•■/  UP  GENERATES  AM  INPUT  CORRESPOND I k'G  TO  A 

e2  C*  ■••■***  (UNIT)  STEP  I"  VELOCITY.  THIS  IS  A  (UNI  T)  RAMP  IK 

o4  C- ***>  RANGE. 

cb  SUBROUTINE  KAMP(Z.T) 

cc  uI.ME'.SKi:  Z(2) 

c7  Z( I )=Z( l)+Z(2)*T 

cu  Z  (  2 )  =  I  .  1. 

CV  Z(3)=U. 

7 1:  licVJii.T 

7 1  e:ii: 

^  *•.>*«■  *  tVix*  •A-'*  *>  %A  A  ArA  *  ArAr * ArAr  ArAr A  *  ?VAr  ★ 

vs  c-' i-..s:;;<K  yjtine  pa  rad  generates  am  input  corresponding  to  a 

7<  b* ( ij { : I  T)  ilEP  I!!  ACCELERATION.  THIS  IS  A  (UNIT)  PARABOLA 

7  b  cv*.<*T.-r;  range. 

7  o  SUBROUTINE  P/N  A.F3  (Z  ,T  ) 

DIMENSION  ZC) 

7  c,  Z(  I  )=Z(1  )+Z(2)*T+Z( 3  )*T**2 

7  y  Z(2)=Z(2 )+Z(3)*T 

2..  Z(3)  =  l.i. 

I  RETIJNN 

.;2  Em; 

(Jt  A.  ;V *.!; A  AAA.  A.  A  >  A  >  A  A  A  A ★★  A:  Ar  ★★  A  *  *  Ar  Ar 

2  4  c*^**^SUL-KOUTi;:c  ZERO  INITIALIZE?  parameters  prior  to  ITERATING 
Ob  t;*p***» THE  ICAL-AI!  ALGORITHM  BY  SETTING  THE  STATE  VECTOR  (Z) 
a  i>«.**»A.NP  THE  PREDICTOR  VECTOR  (Zl'  I G)  TO  ZERO.  IT  ALSO  ZEROS 

1.7  disc,  parameters  necessary  to  the  computation. 

I  -•  SUBROUTINE  Z  Elt  0  ( Z ,  Z'.'i  I G , ER  !  A  X ) 

c >  DIMENSION  Z(2)  ,Zt.TG(3) 

V.  NO  I  v-  1  =  1,3 

VI  Z(I)=0. 

V2  Z'..IG(I)=i. 

V2  10  CON  ITUit 

C  Hr.,'  A  X=!l. 

v  b  2.  E  i  Ur!'. 

VC.  tNI' 

V  7  0  V  V  *??>,  V  vc  *  ■:  *  *  *V.  V  *******  ******  *** 

v  r  i_  •'.***•..  SuHhOUT  I NE  ERROR  1  COMPUTES  THE  MAXIMUM  ERROR  tOVFRSHOOT) 
VV  C  *  *  *■•  r  UR  tf  A.1PS  AND  PARABOLA  INPUTS 
li.  SU  ’h'JUTINE  ERROR  I ( Z.ZUIG, ERR  ,  ERI’AX  ) 

I!  i  i:RR=Z-ZI  IG 

i,..-:  a;jer,c=ai.S(ERR) 

IN-  1 1- (El;  IAX.LT. A3ERR)  ERMAX=APERR 

lv  *  RllTUR:! 

1 1  b  END 

i>  C  o*  **  *“*  v  *****  *s.  ****************  * 

1.7  0* **** SUpR  OUTI NE  ERR0R2  COMPUTES  THE  ABSOLUTE  VALUE  OF  THE 
1.  0  C*»**»KFIo'n  OF  ‘THE  FIRST  OVERSHOOT  IN  RESPONSE  TO  STEP 

1 1  v  0*  ****CI  IANt.cS  I  N  RANGE. 

II;.  SUBROUTINE  E  RR0R2  ( Z ,  ZV.T  G ,  ERR  ,  ERM  AX  ) 
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I  i.)  <-.*  Vf'.V*,V  v  X  ..  r'f5.  ¥:**■*•  7  *  x  *  V  v.*  *** 


i  i  v  i.****»sy ,.  cVi **::  ok  kh'tes  11::  oof  step  kaljov:  prediction 

ij.  here  thi:  variable  is  the 

i.:  i  v.**-.  ,\ r  ■  ::',a  1,3  m:  i"\'ovatio?:. 

i  sk::  :  ;<3  >’op<3> 

I.:-  r::  ; 

...i  (-•> 

i>  ( i ; 

us  ...i  < :  )=;:.!  '( i  )+'■  r'<2>*T(  n+z;  lo<3>*T<2>/?.*'3r(  i  )*ir.ov 

13-  let.')-:..  I.  <3  )<•!;.  10(3  )*T(  1  )*iJ»(2)*r'-:0V 

i-  •.  ..  ..  I. A.>  )  =  k  (3)+w"(.3)<rir:  o\ 


..  .  /jvi:  ,v.::;s  t -<a: sf:-:i;s  the  pairs  feci  the  array 

i : ■'  v-. ..oa 1 1.  Tt.  j he  v kctck  op 
>:■,  !«;3Ti:  •:  i :-’h < xha i ?: . •  ;r> ) 

I--  1  :  :e..sic.‘:  xo/p'O) ,gp< : ) 

i-'.  i  ;  =  .<(/;•.(  t ) 

i:;  c.-?(2  )=,<i  aii'(2) 

i3*  jp(3)=:;c/i;:(3) 

i*..  ret:::.:; 

i  *  t  .' 

I  i  2  L*  *<*■■.  *:»******  **  **********  ***** 

i-_  !_»**»*:■  o'  ,3Ti:i3  STOP  CHECKS  TO  SHE  IF  THE  YAXIfW  ERROR  H AS 

i * -•  c*.  •  :. he" :  Aci-.u-vnn.  if  r-c,  the-  a  flag,  istcp,  is  set. 

t-v  ?3brtOCTl :=  STOP(E!(E.  EK'IAX,  ISTOP) 

1  *■'  c  I  STOP=i ! 

:**.  AbE«::  =  AFS<  EUR ) 

!•■  If  CASb’lil.  .I.T.  EIH'AX)  ISTO°=  i 

'  -  *  ru:  I  OR!! 

i  o .  r. !  f  i . 

I  0  i  0*  **-..*.  V  ***  :r**>  **********  ******* 


Appendix  VI 

The  Computer  Program  TRANS 


1  o*  **** **** *★  P  ROGK A' *  TRANS  ****★*★★★***•.?*★ 

2  opt  i  o:  is  op 
include  •inssF.SYsv 

4  DlKeiJSIlK  Z(2>  ,Z'.iIG(2)  .GP(2)  .P/IG(6)  1 00 ( 6  ) 

5  Di/.'.ENSicr  xav(<i .4 1 .5) .yoijtu, <1  .5) 

C  HEAL  LA/,  t  0  AC  6  > 

7  COMIOM/Fill/RHO.l  (4  ) 

t  COM/! 01 J/SET2/P I »V .  V  All  V ,  S I Oil ,  V  ARJ  .  SIGA.7ARA 

0  OAT  A  L  ALL  DA/  I . ,  .  5,  .  1  , .  05 , .  0 1  / 

I  u  CALL  ASSIGII(6HTKN0C7,‘jU4U)3S.I  ) 

11  CALL  ASSIGN  ( <5MTCT0( ‘7  ,5r!4  I03S.2  ) 

12  CALL  DEASSN 

13  CALL  SUPEkR<3) 

I  4  T  <  I  )  = .  0  I 

IP  T(2 ) =1 ( I )**2 

1  o  T  ( 3 )  =T  (  I  )  **3 

17  T( 4 ) =T( I ) **4 

Id  S'  I GA  =2  • 

IV  V ARA=5I CA**2 

2C  r  ACTCu=  10.  **  .  I 

2  1  CO  950  ICfiTfiL=  1 .2 

22  co  y.io  i  i.a/’=  i ,  ;i 

22  FliiO=l  .-LA'.IBD/(  ILAM)*T(  I) 

24  k*tt*FAC  i  CR 

2-j  R=U;.**-5/FAC70f( 

2c  DO  HUM  IPX*  1, 4  I 

27  K=K*FACTOH 

2b  SI  CV  =R*S  IGA/LA /7!DA(  I  LAM)  **2 

2  V  VARY  =SI CV»*2 

21.  VAH.I=(  I  . -rlHO**<?)*V Al/A 

2  1  SIGi*=SOkT(VAf<:.'> 

22  CALL  ZERO ( Z, ZN I G,  ERNAX.K ,  I  5  TCP 

22  CALL  SET'1P2(  P-ilG  ) 

24  DO  7 Otl  1  =  1  ,5fc'v)C) 

25  CALL  GPCALC(PRIG.CP) 

2 c  GO  70  (  U1.2C.30)  .ICNTRL 

-7  10  CALL  STEP(Z) 

2b  CALL  EKR0R2 ( Z(  I  )  .Z'AIOC  I  )  ,ER'H  .ERMAa  ) 

2  V  GO  TO  5C 

40  20  CALL  RA/.;P(X,T(  I  )  ) 

4  1  CALL  ERROR  I ( Z(  I  )  ,ZV!IG(  I  )  .ERR ,  Efi*' AX  ) 

42  GO  3'0  5C 

42  20  CALL  PARAB(Z,T(  I  )> 

‘4  40  CALL  ERROR  I  ( Z(  I  )  , ZUI  0(  I  )  ,  ERR, ERNAX ) 

45  SO  CALL  COUNT ( I, K, ERR) 

40  CALL  STCP2  (Z.ZI'IG.I STOP ) 

47  CALL  PRED(Z(  I  )  .Z'-'IG.GP) 

4b  CALL  PUKPWIG.PV.TGO) 

4y  IrdSTOP.EO.  I  )G0  TO  ICO 

50  7 CL  CONTINUE 

51  ICC  XOUT ( I , I PT, I  LAM) =H 


52 

52 

5«  oC 
55 


YOU! ( I  ,  IPT.ILA  l)=R 
GO  TO  (OM.70.HO) .ICUTRL 
XOUT (2,1  Pi , I  LAV ) =EKNA  X/. 20 
Y OUT  ( 2 , 1  PT ,  I L  A '  I )  =K 
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Z  ..iu<  I  )=>■). 
ii.  CONTINUE 
ERMAX=i'>. 

K=D 

I  SI  Cf>=.- 
RETURN 
Eillx 

(_>  *:<3r»  *****>ii 

C**»*»SUU.<OUVl NE  EKiiOli  I  COMPOTES  THE  .MAXIMII'  ERROR  (OVERSHOOT) 
c*****rOR  ramps  Ann  parabola  inpUxS 

SUBROuYl  !.’E  E  RROR  t  (  7. .  A' 1 0 ,  ERR  .  ERM  A X  ) 

EnR=Z-Zx IG 
ABETn<=Al  S(EtvR) 

IKfc'H’IAX.LT./NEEK)  EKHAX-ATERt 
mExT'ki'J 
EMD 

XXXX  xxxx  XX  XX  xxxx  xxxx  xxxxx 

l*xxx-.sugkouvi::e  Ekrorz  ec.  putts  ti«e  absolute  value  of  the 
c****»i!Eio:n  u-  tub.  first  overshoot  it  res^nse  to  step 

C » x xx* CH Ai (E S  III  H A’.'Oli . 

SU BE  CL1  x  I  f  'E  E  RECrx2  (  7 «  ZN I  G  *  b  i<T? ,  El./1  A X  ) 

EEH=Z-ZMG 

I r  ( fcmt . L'x .  C .  )OC  TO  1  : 
t(EToE!J 

it  AUEr.K=AL  S(  ELI: ) 

I  r  ( b'RMA.V.L  X .  A  DERR  )  ERF  AX=  A  BEER 
KhTUrXN 

END 

lie*  k *  ★  *  **  x-*  **  *  ★**■*  *  Hr  sir  ir  k*'*  *  *  *  *  * 

C***x* SUBK0UT1 1;=  PEED  COMPUTES  Y’-'c  ONE  STEP  KAIT'AN  PREDICTION 
txxxx»liS  IMG  ECU  ATI  ON  3.  HERE  TIE:  VARIABLE  X  IS  THE 
C*  xxx*  I IIPU7  AID  It-TOV  is  THE  l  CMC  VATIC”.. 

SUBROLi  XT  NE  PE.EDIX.Z:  IG,GP> 

DIMENSION  Zl.TGO  ) .  GP  (3  ) 

REAL  I  Nt.CV 
CCi.f.iOU/VT'1  /RFO  *T(4 ) 

INNOV=X-ZN!C( I  ) 

ZS/'IG (  I  )-Z  .TG(  I  )  +  Z'.TG(2  >*T(  I  >+Z.  IG(  3)*'.  <2  )/.?  *+GP(  I  >*I”.’  OV 
ZU IG(2 )=Zn IG(? )+ZN 10(3 )*T(  1  )  +GI  <  2 )  *1 NUOV 
Z  l .  I G  ( 3  )  =  k  H  0*  Z  id  G  ( 3 )  +  G  P  (  3  )  *  I ! '  HO  V 
WETUmN 
END 

|  [J*  ***%.********  ******  ****  ****  **  * 

IS4  C» **XxSN-BRO'JTKiE  GAIN’S  TRANSFERS  THE  GAINS  FRO’.'  THE  ARRAY 
ItS  CxxxxxaGAIN-  IG  THE  VECTOR  GP 
ISo  SUBROUTINE  GA I  MS  C  XG  A I II  ,GP) 

It'.  Ill  .(ENSIGN  XGAI  N't  3 ) ,  GP(  3 ) 

\-y.j  GP  (  1  >=XCA  1 1.(1) 

itv  GP ( 2 ) =XGA I N ( 2 ) 

ic.  Gi’(3 )  =xg  a  i  i  :  ( :■ ) 

IC-I  RETURN 
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lob 
12c 

I  4  3 
)  40 

13V 
I  40 
14  I 
142 
1 4  3 
1  44 
I  4b 
l4o 
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4*  xxxx  :*)2t:))  A  Cl  IcVEt'*  Ir  ?c*  . 
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IF  THE  MAX  I  JUKI  ERROR 
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22  I  cut  or ! /  i  H  i  2 /  H  ( )  V  ,  V  A M  V ,  S I  ;  ,  V  A W IV ,  S I  n A  .  /  A  l(  A 

2.2  i-.i  (  I  )=5.*i/Ai(V+  {  |  .  +  :;H(J)*VAH/*T(4  )/2  . 

22-  PA  (2  )  =  J.*(  VAliV/Yf  1  )+(  I  .  +  2HC  )  *v'  AK  A*  1’  (  3 )  /+  . ) 

22  P  m  (3>=(  1  .+m!C)*liH0*VAHA*T(2  ) /?. . 

2+0  (  A  )=2  .  *V  APV/T ( 2  )  +(  1.2+  <42 )  W  Arf  A*T(  2  ) 

2+o  Pm  ( 0  )  =  (  I  . +HHC/2 .  )* 2 Al( A* T  (  I  ) 

2+  1  (’1.  ( o  )=s<HC'**2*V  Aif  A 

22  »  tifl  UitN 

22  V  LM1 

2  2  2  ^xXxxXXXXXXXX*  X  X  :V  X  X  XXX  XX  XX  X  X  XX  X 
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22+  ul  'JblU:  PUol.GPC  ) 

22+  C0-M'.':l/l-rI/iir0,T(4) 
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2+o  p=  I  . /(  PI.  (  I  )  + V  ARV  ) 

2+1  uP(  I  )  =  (!'.  (  I  )  +  PI-  <2)-n  <  I  )  +  !-.•.  (+  ) *  1  (  2  )  /2  . )  * 

2  3B  GP(2  )  =  (Ff  (2)+Pl.(3)*T(  I  )>*P 

2+s.  GP(3>=PI  ( 3  ) * (il •  0* P 

2  +  2  KLTUrfiJ 

2+1  tHI.' 
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